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vant Schwarzschild-Milne equation is solved exactly in the absence of internal reflections. 
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1. INTRODUCTION 



The theory of multiple light scattering has been a classical subject of interest for 
one century, which attracted the attention of many scientists, including Lord Rayleigh, 
Schwarzschild, and Chandrasekhar. Standard books are available, such as those by Chan- 
drasckhar [1], Ishimaru [2], van de Hulst [3], and Sobolev [4]. The discovery of weak- 
localisation effects, and chiefly the enhanced backscattering cone [5], yielded a revival of 
theoretical and experimental work in the area of multiple scattering in disordered media. 
Much progress has been done recently in the analysis of speckle fluctuations [6] , in analogy 
with the conductance fluctuations observed in mesoscopic electronic systems [7]. 

Laboratory experiments are often performed either on solid Ti02 (white paint) sam- 
ples, or on suspensions of polystyrene spheres or of Ti02 grains in fluids. In most cases, 
the wavelength Aq of light in the diffusive medium, the scattering mean free path £, and 
the thickness L of the sample obey the inequalities Xq <^ £ <^ L. Multiple scattering 
is also of interest in biophysics and medical physics, in order to understand the trans- 
port of radiation through human and animal tissues. Besides light, all kinds of classical 
waves undergo multiple scattering in media with a high enough level of disorder, i.e., of 
inhomogeneity. Well-known examples are acoustic and seismic waves. The propagation 
of electrons in disordered solids also pertains to this area, since Quantum Mechanics also 
basically consists in wave propagation. 

Multiple scattering of waves in disordered media admits the following three levels of 
theoretical description: 

(i) The macroscopic approach consists in an effective diffusion equation, which describes 
the transport of the diffuse (incoherent) intensity I(r, t) at point r at time t. This approx- 
imation turns out to be very accurate in the bulk of a turbid medium, and more generally 
on length scales much larger than the mean free path i. The diffusion approximation 
yields several interesting predictions, among which we mention the 1/L-decay of the total 
transmission through an optically thick slab of thickness L ^ £, the decay time of the 
transient response to an incident light pulse, or the memory of a typical speckle pattern 
when the frequency of light is varied. This approximation also allows for a quantitative 
prediction of the diffuse image of a small object in transmission [8] . 

(ii) The mesoscopic approach, used by astrophysicists throughout the classical era of the 
subject, is referred to as radiative transfer theory [1-4]. This theory relies on the radiative 
transfer equation, which is a local balance equation, similar to the Boltzmann equation in 
kinetic theory, for the diffuse intensity 7(r, n, t), with n being the direction of propagation. 
This approach leads in a natural way to distinguish between the scattering mean free path 
i and the transport mean free path £*, to be defined below. The diffusion approach (i) is 
recovered in the limit of length scales much larger than i*. 

(iii) The microscopic approach consists in expanding the solution of the wave equation in 
the disordered medium in the form of a diagrammatic Born series. In the regime i ^ Xq 
the leading diagrams can be identified, in analogy with e.g. the theory of disordered super- 
conductors [9] . For the diffuse intensity they are the ladder diagrams, which are built up by 
pairing one retarded and one advanced propagators following the same path through the 
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disordered sample, i.e., the same ordered sequence of scattering events. This picture agrees 
with that behind the radiative transfer equation, which is a classical transport equation 
for the intensity. The ladder diagrams can be summed and yield an integral equation of 
the Bethe-Salpeter type for the diffuse intensity, which we refer to as the Schwarzschild- 
Milne equation. The radiative transfer approach (ii) is thus recovered in the weak-disorder 
{£ ^ Ao) regime. A further step consists in including some of the subleading diagrams, 
of relative order Xo/£ ^ 1, which account for interference effects between diffusive paths. 
Among those contributions, the class of maximally-crossed diagrams is of particular in- 
terest, since it describes the aforementioned enhanced backscattering phenomenon [10, 
11]. 

To summarise this discussion, each of the descriptions mentioned above represents a 
qualitative improvement with respect to the previous ones (see e.g. ref. [12]). In order to 
derive quantitative estimates of observables in the regime £ » Aq of interest, it is sufficient 
to consider radiative transfer theory. The macroscopic approach (i) is clearly insufficient, 
since we aim, among other things, at a description of the crossover from free propagation 
to diffusive transport which takes place in the skin layers of thickness of order £, near the 
boundaries of the disordered medium. This phenomenon is, by its very definition, beyond 
the scope of the diffusion approach. The radiative transfer approach (ii) leads to study 
the Schwarzschild- Milne equation (2.26, 28). This equation has only been solved exactly 
in a limited number of cases. Analytical results have been obtained for isotropic scattering 
of scalar waves [1, 13-15], and in the case where the scattering cross-section depends 
linearly on the cosine of the scattering angle [1, 16]. The case of anisotropic scattering has 
essentially been investigated by means of general formalism and by numerical methods [2, 
3, 17]. 

For a finite sample, physical observables such as the reflected or transmitted intensity 
in a given direction depend on the particular realisation of the sample under consideration. 
Such quantities arc indeed the results of intricate interference patterns throughout the sam- 
ple; they only become self-averaging quantities in the limit of a large enough sample. This 
definition can be made precise by means of the dimensionless conductance g ~ N(./L, 
related to the number jV ~ A/X^ of open channels, where A is the transverse area of 
the sample. The self- averaging regime corresponds to ^ ^ 1. The whole distribution of 
observables is therefore of interest, as long as g is not very large. We mention a recent 
experiment [18], where the third cumulant of the total transmission, an effect of relative 
order 1/g^, has been measured and compared to theoretical predictions. This third cumu- 
lant is of the same order of magnitude as the universal conductance fluctuations, either in 
electronic [7] or in optical [6] systems. In fact the full distribution of the total transmission 
through an optically thick slab has been derived recently [19]. In the following we focus 
our attention on the mean values of physical quantities. 

The vector character of electromagnetic waves also introduces its own intricacy. In 
general four coupled integral equations have to be solved, which are associated with the 
Stokes parameters of the diffuse light in the medium. These equations have been solved 
exactly in the case of Rayleigh scattering, i.e., the regime where the size of the scatterers is 
much smaller than the wavelength [1, 20]. Among specific features pertaining to diffusive 
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light propagation, let us mention the dependence of the backscattering enhancement factor 
on the polarisation states of the incoming and outgoing radiation [21-23], or the progressive 
destruction of the backscattering peak induced by a magnetic field, due to the Faraday 
rotation in a magneto-optically active material [24-27]. 

Furthermore, in practical situations the optical index uq of the scattering medium is 
often different from the index ni of the surrounding medium. This index mismatch, mea- 
sured by the ratio m = no/ni, causes reflections at the interfaces. In the regime of a large 
index mismatch (m <^ 1 or m 3> 1), the transmission across the interfaces is very small, 
so that the light is reinjected many times in the diffusive medium. As a consequence, the 
skin layers become very thick in this regime. More generally, the diffusion approximation 
works better and better as the index mismatch gets large. Improvements of the diffusion 
equation have been proposed [28, 29], which take internal reflections into account. It is in 
fact possible to derive analytical expressions for the reflected and transmitted intensities 
and other observables in this regime. This asymptotic analysis has been performed in 
ref. [15] for isotropic scattering. It will be generalised hereafter to the case of general 
anisotropic scattering. This approach provides accurate results, even for a moderate index 
mismatch m. 

More generally, one of the main goals of this paper is to quantify the dependence 
of quantities on the anisotropy of the scattering mechanism. It has been known for long 
from radiative transfer theory that two length scales are involved in the case of general 
anisotropic scattering: the scattering mean free path i, namely the effective distance be- 
tween two successive scattering events, and the transport mean free path £*, namely the 
distance over which radiation looses memory of its direction. Both mean free paths, to be 
deflned more precisely in section 2, depend on details of the scattering mechanism, such as 
the shape, size, and dielectric constant of the scatterers. In the regime of very anisotropic 
scattering, we have i* ^ i. The ratio r* = £* /£ ^ 1 will be referred to as the anisotropy 
parameter. In some experimental situations r* can be of order ten or larger [30]. The 
regime of most interest is then Xq <^ £ <^ £* <^ L. 

The setup of this paper is as follows. In section 2 we present some general formalism 
on radiative transfer theory and we derive the associated Schwarzschild-Milne equation. 
We show how solutions of the latter equation yield predictions for quantities of interest, 
such as the diffuse reflected and transmitted angle-resolved intensities. The determination 
of the shape of the enhanced backscattering cone is also addressed. Section 3 is devoted 
to the regime of a large index mismatch, for general anisotropy. In section 4 we derive 
a complete analytical solution of the radiative transfer problem in the regime of very 
anisotropic scattering (r* S> 1), in the absence of internal reflections. The paper closes up 
with a discussion in section 5. 

2. GENERALITIES ON ANISOTROPIC MULTIPLE SCATTERING 

Throughout the following we restrict the analysis to multiple scattering of scalar waves 
by scatterers located at uncorrelated random positions, in the regime where the scattering 
mean free path £ is much larger than the wavelength Aq of radiation in the medium. We 
have summarised some useful notations and deflnitions in Table 1. 



4 



As stated in section 1, we put special emphasis on the dependence of physical quan- 
tities on the anisotropy of the scattering mechanism. After averaging over the random 
orientations of the individual scatterers, the differential scattering cross-section of arbi- 
trary anisotropic scatterers can be written in the following form [1, 2] 

da{n, n') = {u/47Tfp{Q)dQ' . (2.1) 



In this formula u is the scattering length, n and n' are unit vectors in the incident and 
outgoing directions, O is the angle between these directions, so that cos© = n.n', and 
dfl' is an element of solid angle around the direction n'. We assume that there is neither 
absorption nor inelastic scattering. In other words the albedo is unity [the situation of 
non-conservative scattering will only be considered in section 2.7]. The phase function 
p{Q) then obeys the normalisation condition 

The total cross-section reads a = u^/ (47r), and the scattering mean free path £ is given by 

<=^ = i^, (2.3) 

where the density n of scatterers is assumed to be small, in such a way that we have 
Ao. 

In the particular case of isotropic scattering, the phase function p{Q) = 1 is a constant. 
In the general case of anisotropic scattering, the phase function p{Q) is a non-trivial 
function of the scattering angle. As recalled in section 1, one has to distinguish between 
the scattering mean free path i, which is the typical distance between two successive 
scattering events, and the transport mean free path £*, which represents the distance over 
which radiation looses memory of its direction. Both mean free paths are related by the 
following expression, well known from kinetic theory, 

^*^J^ i-(cose)' ^'^'^^ 

with ^ 

(cosG)= / ^^^cosep(e). (2.5) 



The dimensionless ratio r* will be referred to as the anisotropy parameter. We have 
usually T* > 1, i.e., £* > £, since the scattering cross-section is often peaked in the forward 
direction, in a more or less prominent way. The regime of very anisotropic scattering, 
where p{Q) is strongly peaked around the forward direction, corresponds to £* ^ £. This 
regime will be investigated in detail in section 4. 
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2.1. General formalism 

In this section we present some general formalism on anisotropic multiple scattering, 
extending thus the treatment of the isotropic case presented in ref. [15]. Some of the 
results exposed below are already present in lecture notes by one of us [12]. 

We begin with a reminder of radiative transfer theory. In the regime £ ^ Aq under 
consideration, in the absence of internal sources of radiation, and in stationary conditions, 
the quantity of interest is the specific intensity /(r, n) of radiation at the position r, 
propagating in the direction n. The specific intensity obeys the time-independent radiative 
transfer equation, that takes the following local form 

£n.V/(r,n) = r(r,n) - /(r,n). (2.6) 

The quantity 

r(r, n) = y ^p(n, n')/(r, n') (2.7) 

is commonly referred to as the source function. We recall that the phase function p{n, n') = 
p{Q) only depends on the scattering angle 0. 

As recalled in section 1, the radiative transfer equation (2.6) [1-4] can be considered 
as a mesoscopic balance equation for the light intensity inside the diffusive medium, some- 
what analogous to the Boltzmann equation in the kinetic theory of gases. It is equivalent 
to the Bethe-Salpeter equation, obtained by summing the ladder diagrams of the Born 
series expansion of the intensity Green's function of the problem. These diagrams are the 
dominant ones for £ ^ Aq, i.e., to leading order in the density n of scatterers. 

We consider a sample of diffusive medium in the form of a slab of thickness L, limited 
by the two parallel planes z = and z = L. The mean optical index no of the slab can be 
different from the index ni of the surrounding medium. The index mismatch, measured 
by the ratio m = no/ni, generates internal reflections at the interfaces. We introduce the 
optical depth t = z/£ oi a, point in the sample, and the optical thickness h — L/i oi the 
sample. Finally we use angular co-ordinates as in Table 1: ^ is the incidence angle, with 
the notation = cos 9, while the azimuthal angle is denoted by ip. 

Since the problem has rotational symmetry with respect to the 2;-axis, normal to 
the sample, and translational symmetry in the (x, j/)-plane, it is natural to express the 
(^-dependence of the specific intensity and of the source function as Fourier series of the 
form 

/(r,n)= I^"'\r,f^)e'^^, r(r,n)= J] T(^\t, f,)e'^'^ , (2.8) 

— oo<m<+c» — oo<m<+oo 

where the integer m is the azimuthal number. 

Furthermore, along the lines of refs. [1, 4], for general anisotropic conservative scat- 
tering, we expand the phase function in Legendre polynomials as 

p(e) = J2^2£ + l)wePi{cos 6). (2.9) 
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We have wq = 1 (see below), while the other coefficients are only constrained by the 
positivity of the phase function. 

In co-ordinates related to the sample, the phase function then reads 

p{n, n') = if, /x', (^0 = Yl Prnil^^ /x')e^"^(^-^'), (2.10) 

— oo<m<+oo 

with 

Pmil^,!^') = Yl (2^+ nl}nm(/^)^^,rn(/x')- (2.11) 

£>\m\ ^ 

The radiative transfer equation (2.6) thus reduces to 

^i_/M(^, /.) = r(-)(T, /.) - /('-)(t, /.), (2.12) 

or equivalently 

At-^[/("^)(^,At)e^^1 =r("^)(T,/x)e^/^, (2.13) 
and the source functions T^'^\t, ji) are related to the specific intensities I^'^\t, ji) through 

rM(T,//) = p(^)[/("^)(T,//)], (2.14) 

where we have introduced the integral operators 

©(-)[$(/.)] = ^Pn.(/x,/x')$(/xO- (2.15) 

We mention for further reference that the Legendre functions {P^,rn(/^), ^ > 1^1} form 
a complete set of orthogonal eigenfunctions of the integral operator 'D^'^\ with eigenvalues 
vDi, i.e., 

p(-) [p,,M] = mP£,M. (2.16) 

Especially for m = we have Pi,o{id) — -P^(/x), where the Legendre polynomials P(,{jJi) 
already appeared in eq. (2.9). 

The following Legendre functions will be of special interest hereafter 

PoM = PM = l, PiM = Piif^) = f^, Pi,i(//) = yr^. (2.17) 

Indeed the identity (2.16) has the following two special cases of interest 

• The eigenvalue wq = 1, associated with the constant eigenfunction Po(a*) = 1 of V^^\ 
is a consequence of the conservative nature of the scattering mechanism, yielding diffusive 
behavior in the long-distance limit. More explicitly. 




(2.18) 
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• The first non-trivial eigenvalue wi of both operators V^^^ and T>^^^ is related to the 
anisotropy parameter r* of eq. (2.4). Indeed, since -Pi(//) = we have 

—po (//, iJ,')iJ,' = wifi, (2.19) 

hence 

wi = (cos e) = 1 - T* = Y~^- (2.20) 

We also notice that all the eigenvalues of the operators D^"^) are trivial in the particular 
case of isotropic scattering, since wq = 1, while tu^ = for £ > 1. 

2.2. Schwarzschild-Milne equation 

We now turn to the derivation of the Schwarzschild-Milne equation in the general 
situation of anisotropic scattering. This key equation of radiative transfer theory will be 
the starting point of the following developments. 

It turns out that most observables of interest can be derived by considering quantities 
with cylindrical symmetry around the ^-axis, namely those corresponding to an azimuthal 
number m — 0. Henceforth we restrict the analysis to this sector, except in section 2.6, 
and we drop the superscript (0) for simplicity. 

We consider first a half-space geometry (6 = oo). We assume that the limiting plane 
T = of the sample is subjected to a cylindrically symmetric incident beam, characterised 
by an angle of incidence 9a, i.e., that the incident intensity does not depend on the az- 
imuthal angle This is automatically satisfied at normal incidence {9a = 0). Under 
these circumstances, the inward intensity on the limiting plane r = 0+ contains both the 
normalised refracted incident beam, coming in a direction defined by /x^, and the intensity 
coming from the bulk of the medium, after being reflected once at the interface. Using the 
radiative transfer equation (2.12, 13) for m = 0, we thus get 

/(0+, fia) = 25ifi -^ia) + ^ r dTe-^/'^r(T, fia) (/^ > 0), (2.21) 

Jo 

where the Fresnel reflection coefficient R{fi) is given in Table 1. The only contribution to 
the outward intensity on this plane comes from the light rays which have experienced at 
least one scattering event, namely 

1 f°° 

/(0+, -II, lia)^- dTe-/'^r(T, -//, lla) (/^ > 0). (2.22) 

^^ Jo 

The radiative transfer equation (2.12, 13), together with the boundary conditions 
(2.21, 22) at r = 0"*", can then be recast in the integral form 

/(r, lia) = 25{fi - //a)e-"/''° + {K * T){t, //«)• (2.23) 
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Here and throughout the following, the star denotes the convolution product 



{K * r)(T, /X, fia) = dr' j ^ ^K{t, r', /x')r(T', (2.24) 

The kernel K can be split into two components: K = JCs + Kl. The bulk kernel Kb 
contains the contributions to the intensity at depth r arising from the scattering from 
either smaller or larger depths. The layer kernel Kl takes into account the intensity being 
scattered at depth r' in the direction of the wall, then reflected there, and then scattered 
at depth r. These kernels read explicitly 

Kb{t, f,, t', = 26{f, - f,')9{fx{T - rO)^e-(---')/^ 

R(a) ^^-^'^ 

The final step consists in using eq. (2.14) in order to derive from eq. (2.23) the 
following closed integral equation for the source function 

r(r, //, ^a) = Po(/^, //a)e-"/'^» + (M * r)(r, //, //„), (2.26) 

which we refer to as the Schwarzschild-Milne equation of the problem. 

The kernel M has the following two components, in analogy with the kernel K from 
which it derives 



(2.27) 



I 

M,(r,/.,r',/.') = ^( V)^^^^^^i?( V)e(^+^')/'^', 

lA* I 

so that eq. (2.26) reads explicitly 

+ r^^' /'^Po(/x,/x')e-(---')/'^'r(r', /.',/.„) 
Jo Jo ^/^ 

+ dr' ^poili, -/.')e-(-'--)/^'r(r', V, /^a) 
Jo Jo ^f^ 

In the case of isotropic scattering, the phase function po{ij,, n') 1 is a constant, and 
the source function T{t, /i, /la) does not depend on The Schwarzschild-Milne equation 
(2.26, 28) thus takes a simpler form, that has been extensively studied in ref. [15]. The 



(2.28) 
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rest of this section is devoted to an extension of the results derived there to the general 
case of anisotropic scattering. 



2.3. Solutions of Schwarzschild-Milne equation and sum rules 

In the general case of conservative anisotropic scattering, the Schwarzschild-Milne 
equation (2.26, 28) has a special solution Ts{t, /j,, ^a) which remains bounded as r — > oo, 
whereas the associated homogeneous equation has a linearly growing solution TniTjH). 
More precisely, it can be checked by means of eq. (2.18, 19) that both source functions 
and the associated specific intensities have the following asymptotic behavior for large 
depths (r 1), up to exponentially small corrections 



The quantities tq and Ti(//a) are unknown so far. Let us anticipate that they play a 
central role in the following, in the sense that they will bear the full non-trivial dependence 
of quantities on the scattering mechanism. These quantities also obey two groups of sum 
rules, (2.37, 38) and (2.41, 42), to be derived below. 

To do so, it is most convenient to introduce the Green's function Gs(t, fx, r', fx') of the 
problem, along the lines of ref. [15]. It is defined as the solution which remains bounded 
as T — >^ oo of the equation 



which merely express the time-reversal symmetry of any sequence of scattering events. 

As a consequence of eq. (2.26), the special solution T s{j, Ha) can be expressed in 
terms of the Green's function as 



Ti? (r, iJ,)^T- h{t* - 1) + Tor*, 
Inir, fJ^^T- fXT* + tqt*. 



(2.29) 



Gsir, fi, t', fx') = poifx, fx')6{T - r') + (M * Gs){t, fi, r', fx'). 



(2.30) 



The kernel K and the Green's function Gs possess the symmetry properties 



Gs{r, II, t', ix') = Gs{t', r, -//), 



(2.31) 




(2.32) 



We also define for further reference the following bistatic coefficient 




(2.33) 
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The latter expression defines the bistatic coefficient for any complex values of its arguments 
with Re > 0, Ke fii, > 0, even outside the physical range ^a-, ^^b < 1- The symmetry 
lilJ'ailJ'b) = lifJ'bj IJ'a) IS a conscqucnce of the properties (2.31). It is thus again due to 
time-reversal symmetry. 

On the other hand, a relationship between both solutions FniTjH) and rs{T, n, Ha) 
of the Schwarzschild-Milne equation can be derived as follows. The Green's function 
Gs{t, fi, t', jj!) is clearly asymptotically proportional to the homogeneous solution F/j (r, jj) 
when r' goes to infinity, namely 

hm Gs(t, /X, r', ii') = ^Th{t, fi), (2.34) 

where the proportionality constant D will be shown in a while to be equal to the dimen- 
sionless diffusion coefficient (2.39). 

As a consequence of eqs. (2.29, 32-34), we have 

T,{fia)= lim 2(^W^= 1 rdTe--/>^'^TH{T,-fia). (2.35) 
^lb^oo /ifo U Jq 

We now turn to the actual derivation of the sum rules obeyed by the quantities de- 
fined so far, which are related to the so-called F and K-integrals, with the notations of 
Chandrasekhar [1]. 

The first group of two sum rules is a consequence of the conservation of the fiux in 
the 2;-direction in a non-absorbing medium, given by the following F-integral 

F(r) = l'^^(,IiT,fi). (2.36) 

It can indeed be checked, using eq. (2.12), that dF/dr = 0. 

We investigate first the F-integral Fs associated with the special solution Ts{t, /U, /Ua). 
Considering the r — > oo limit determines Fs = 0, whereas considering the r ^ limit 
yields the sum rule 

^T{fx)-f{fX,lla) = f^a, (2.37) 

where the transmission coefficient T(//) = 1 — -R(^) is given in Table 1. The jia ^ limit 
of eq. (2.37), together with eq. (2.35), yields another sum rule, namely 

f^^nn)T,{n) = l. (2.38) 

We can evaluate in a similar way the F-integral Fh associated with the homogeneous 
solution Th ij, n) . This yields no independent sum rule, but leads to the identification of 
the constant D of eq. (2.30) with the dimensionless diffusion coefficient 

D = Y- (2-39) 



11 



The diffusion coeflBcient indeed reads -Dphys = c£D = c£*/3 in pliysical units [1-4]. Tfie 
transport velocity indeed coincides witfi tfie velocity of light in vacuum c, to leading order 
in the regime £ ^ Xq. 

Besides the sum rules (2.37, 38), which were already given in ref. [15] for isotropic 
scattering, the radiative transfer equation also admits another group of two sum rules, 
which are novel in this context, and whose intuitive interpretation is less evident. Consider 
the so-called K-integral, again with the notations of Chandrasekhar [1] 

K{r)^ l'^^l,'l{r,f^). (2.40) 

It can be checked that eqs. (2.12, 18-20) yield dK/dr = -F/t*, hence K{t) = -Ft/t* + 
Kq, with Kq being independent of r. By considering the i^T-integrals associated with the 
special solution Ts{t, /U, fXa) and with the homogeneous solution Fj^ (r, //), we obtain after 
some algebra the following sum rules 

+ ^(/^))/^7(/^, /^a) = ^ - (2.41) 
^'^{l + R{^^))^^TM=To. (2.42) 



2.4. Diffuse reflected intensity 

The evaluation of the angle-resolved diffuse reflected intensity by means of the general 
formalism exposed above closely follows the lines of refs. [4, 15]. We consider a half- 
space geometry, and we assume that the limiting plane of the sample is subjected to a 
cylindrically symmetric incident beam, characterised by an angle of incidence Oa- This 
technical assumption includes the case of a plane wave at normal incidence {9a = 0). 
Under these circumstances, the diffuse reflected intensity per solid angle d^lf, reads 

dR{a^b) ^ cosOaTaTi, 

— = A«(6'a, 91,) = 2 7(/ia, (2.43) 

where we have again used the notations of Table 1: m — no/ni is the index mismatch, and 
Ta = T^Ha) and T5 = T(/Xb) are the transmission coefficients in the incident and outgoing 
directions, respectively. 

The essential factor in the result (2.43) is the bistatic coefficient ^{jia,l^b)i whose 
definition and general properties have been exposed in section 2.3. It will be evaluated 
more explicitly in the large index mismatch regime in section 3, and in the very anisotropic 
regime in section 4. 

2.5. Diffuse transmitted intensity 

In this section we consider the angle-resolved mean transmission of an optically thick 
slab, of thickness L = with 6^1 being large but finite. A generalisation of the 
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reasoning of section 2.1 allows us to write down the following Schwarzschild- Milne equation 
in this geometry 

rb(r,/x,/Xa) =po(/^,A*a)e"'^/''" 
Jo Jo 

+ J'dr' ^po{^^, -^^')e-^-'-^y^'n{r', ^2.44) 

+ I'dr' ['^R{^i')po{^^,^i')e-^^+^'^/^'n{T',-^i',^ia) 

Jo Jo 

Jo Jo ^/^ 

The solution of this equation for a thick slab can be constructed from both solutions 
Ts and Th of the half-space geometry, by means of a matching procedure, along the lines 
of refs. [4, 15]. Using the asymptotic forms (2.29), we obtain 

[ Ts{t, I,, i,a) - J'^^"\ '^h{t, y) (r finite, 6 - r » 1), 
r,(r,/.,/.«)^<^ (2.45) 

' 6 + 2^-0!-- ^^^^ " ^' "^^ " ^ ^ ^ 

Both expressions lead to a linear (diffusive) behavior [3, 15] in the bulk of the sample 
(r » 1, 6 — T ^ 1), namely 

The derivation of the diffuse transmitted intensity per solid angle element dftf, again 
closely follows the lines of refs. [4, 15]. We obtain 

with 

A'^{ea,e,) = :^^ri(/.,)ri(/.,), (2.48) 

where we have again used the notations of Table 1: m — no/ni is the index mismatch, and 
Ta = T{iia) and T5 = T^Hb) are the transmission coefficients in the incident and outgoing 
directions, respectively. 

The essential ingredient in the above result is the function ri(/i), whose definition 
and general properties have been exposed in section 2.3. It will also be evaluated more 
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explicitly in the large index mismatch regime in section 3, and in the very anisotropic 
regime in section 4. 

The result (2.47) shows that the effective thickness of the sample is {b + 2tot*)£ = 
L + 2to£*. In other words, zq = tq£* represents the thickness of a skin layer. This quantity 
is also referred to as the injection depth, or the extrapolation length of the problem. 

2.6. Enhanced backscattering cone 

The general formalism exposed above can be extended to the study of the enhanced 
backscattering phenomenon, which takes place in a narrow cone in the vicinity of the 
exact backscattering direction [5] . This phenomenon is due to the constructive interference 
between any path in the medium and its time-reversed counterpart. One of the goals of 
this section is to derive a quantitative estimate of the width of the cone, with emphasis 
on its dependence on the anisotropy of the scattering mechanism. As recalled in section 1, 
the form of the backscattering cone can be predicted by summing the so-called maximally- 
crossed diagrams [10, 11]. This can be performed by means of a careful treatment of 
radiative transfer theory [10, 11, 16, 15, 12]. 

We restrict the analysis to normal incidence {9a = 0), and to the geometry of a 
half-space diffusive medium. We introduce the dimensionless transverse wavevector 

Q = q£, (2.49) 

and its magnitude 

Q = qi = koie' = hee > o, (2.50) 

where 9' and 9 are the incidence angles of the outgoing radiation, and ko and ki are its 
wavenumbers, respectively inside and outside the diflFusive medium, according to Table 1. 

Along the lines of refs. [10, 11, 16, 15, 12], the reflected intensity in the vicinity of the 
backscattering direction, i.e., for ^ <^ 1, ki£ » 1, and Q fixed, takes the form 

A'^iQ) ^ ^ [7(1, 1) + 7c(g) - Po{l, -l)/2] , (2.51) 

where 

• The sum of ladder diagrams, 7(1, 1), yields the background reflected intensity in the 
normal direction, in agreement with eq. (2.43); 

• The sum of the maximally-crossed diagrams, 7c (Q), represents the contribution of the 
interference between the sequences of any number (n > 1) of scattering events and their 
time-reversed counterparts; 

• The subtracted third term is the contribution of the single-scattering events (n = 1), 
which are invariant under time inversion, and must not be counted twice. 

We deflne the enhancement factor 

^^^^"^^(0,0)"^+ 7(1,1) • ^^-^^^ 
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It turns out that the peak value of the interference contribution coincides with the back- 
ground term, i.e., 7c(0) = 7(1, 1). Hence the enhancement factor at the top of the cone, 
namely 

nearly equals two, up to the small contribution of single-scattering events. 

We now turn to the actual determination of 7c (Q) [10, 11, 16, 15, 12]. Basi- 
cally, the transverse wavevector Q causes a de-phasing which amounts to replacing the 
pure exponential damping exp(— r///) of unscattered intensity by the complex exponential 
exp (-(l-iQ.n)T///) B ecause of the vector nature of Q, the source functions r("^)(g,T, //) 
pertaining to all sectors defined by the azimuthal integer m are coupled to each other. 

We choose co-ordinates such that Q is oriented along the positive y-axis, in order to 
simplify notations. The source functions then obey coupled Q-dependent Schwarzschild- 
Milne equations of the form 

r(-)(Q,T,/x) = 6m,oPo{f^,l)e-^ + *rW) (Q,T,/x), (2.54) 

— oo<fe<+oo 

generalising eq. (2.26). The Q-dependent Schwarzschild- Milne kernels read M^"^''^) = 
Mb("^''=) + Ml("^'^), with 

Mb (Q, r, r', ^^') = e{^i'{T - r') f-^\'']f\ -^^-^'y^' 

Ia* I 

^ . ' (2.55) 

Mi,''"-''\Q,T,f,,T',f') = B(-^' f t R(-I'')e^^*^'"''' 

1/^ I 



X 



where Jm{z) denotes the Bessel function of integer order. The source functions have the 
following property 

T(-^\Q,T,fi) = (-l)-r(-)(Q,r,/.) (2.56) 
which they inherit from an analogous symmetry property of the Bessel functions, i.e., 

J-n^{z) = i-ir Jm{z), 

Finally, the shape of the backscattering cone is given by 

lc{Q)= dTe-^r(°)(Q,T,-l). (2.57) 
^0 

The top of the backscattering cone, described by the small-Q behavior of 7c(Q)j is of 
special interest, especially because of its universality. Indeed it is due to the contribution of 
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long paths in the diffusive medium, along which the radiation undergoes many scattering 
events. On the contrary, the wings of the cone, corresponding to a large reduced wavcvcctor 
Q, only involve short sequences of two, three, etc. scattering events, and are therefore 
expected to depend on the details of the scattering mechanism. This is already apparent 
in the subtracted term in eqs. (2.52, 53), which involves the single-scattering cross-section 
in the direction of exact backscattering. 

The universal small-Q behavior of the backscattering cone can be determined as fol- 
lows. We look for an approximate solution for Q <^ 1 to the coupled Schwarzschild-Milne 
equations (2.54) as a decaying exponential, with an inverse extinction length sq, times 
expansions in Legendre functions of the form 

r(-)(Q,T,/i) = e-^°"G(^)(Q,//), with G("^)(Q,/^)= 5^ Q,^^A-(/^)- (2-58) 

e>\m\ 

First, we observe that the arguments of the Bessel finictions in the kernels (2.55) are 
proportional to Q. Since we have Jm{z) ~ {z/2)'^ /m\ for small z and m > 0, we therefore 
expect that the coefficients of the expansion (2.58) fall off as C(,^rn ~ Q'"^'. By virtue of 
the symmetry (2.56), we can thus restrict the analysis to the sectors m = and m = 1. 
Second, we make the hypothesis, to be checked later on, that the inverse extinction length 
So is proportional to Q. Then, deep in the bulk of the medium, i.e., for r 3> 1, the integral 
equations (2.54) can be approximated by differential equations, obtained by expanding 
the source functions in powers of (r' — r), keeping only the first two derivatives, and 
consistently the first two powers of Q. We thus obtain the following two coupled equations 



r(0)(Q,r,/.) = p(°) 



F W (Q, r, ^i) - /.Af W (Q, r, ^i) + /x^ /^F ^ (Q, r, ^i) 
dr dr^ 



- ^ (1 - /.2)r(o) (g, r, ^i) + gyr^F(i) (g, r, /.) + 



FW(g,T,/.) = p(^ 



f(^) (g, t, /.) - 1 v'r^F(o) (g, T, /.) + •• • 



(2.59a) 
(2.59b) 



We first solve eq. (2.59b) as follows. Since we only need a leading order estimate of 

F*^^)(g, r, //), we only keep the first coefficient ci^i of the expansion (2.58). Using eqs. 
(2.16, 17), we thus obtain 

ci,i!^-|(t*-1)co,o. (2.60) 

By inserting this last result into eq (2.59a), and making use of eqs. (2.16) and (A. 4), we 
obtain the following recursion relations for the coefficients q^q 

/ £{£-1) 2£^ + 2£-l (£+!)(£ + 2) " ^ ' ^ 
^ \{2£ - 1){2£ - Sf^-^'^ ^ {2£ - 1){2£ + Sf^'"^ ^ (2£ + 3)(2£ + 5) 
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It is clear from the structure of these relations that, when Q is small, the coefficients 
C£,o ~ decay rapidly. Keeping this hierarchy in mind, and using wq = 1, we obtain the 
estimates 

ci,o ~ Q{r* - l)co,o (2.62) 

and 

So ^ Q. (2.63) 

This last result corroborates the hypothesis made in the derivation of eq. (2.59). We shall 
come back to its meaning at the end of section 2.7. 

The next step also follows the lines of ref. [15]. The small-Q behavior of Tc{Q-:T, ii) 
has a term proportional to Q, which is proportional to the homogeneous solution Vh{t^ jj) 
of the Schwarzschild-Milne equation (2.26, 28). Indeed, consider the right-hand side of 
eq. (2.54) for m = 0. The leading Q-dependence there comes either from the action of 
j\^(o,o) Y^^\ or from the action of M*^*^'^^^ on F*^^^-'. All these explicit Q-dependences 
begin with Q^. Putting everything together, we are left with the following estimates of 
the source function r^*^) (Q, r, /x). For Q and fixed r we have 

r(°)(Q,T,//) = Vs{T,^i, 1) - QT^{l)VH{T,^i) + C»(Q'), (2.64) 

whereas for Q <^ 1 and r ^ 1 simultaneously we get 

fW (g, r, /x) = Ti(l)e-Q^ {l + Q(/x(t* - 1) - t^t*) + 0{Q^)^ . 



(2.65) 



The universal peak of the backscattering cone is then evaluated by inserting the esti- 
mate (2.64) into eq. (2.57), using eqs. (2.33, 35). We thus obtain the following expression 

7c(g) = 7(1, 1) (l - ^ + ^(^')) ' (2-66) 



where the width of the cone reads 



i.e., in physical units 



with ki being the wavenumber of radiation outside the diffusive medium. This simple l/£* 
law is already predicted by the diffusion approximation [31, 32]. 

2.7. Extinction and absorption lengths 

Up to now we have assumed that the diffusive medium is conservative. This means 
that the light only experiences elastic collisions; there is neither absorption, nor inelastic 
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scattering, implying the normalisation (2.2) of the phase function. We now want to discuss 

briefly the case of a weakly absorbing diffusive medium, characterised by a non-trivial 
albedo a such that 1 — a <^ 1. In this case the diffuse intensity is expected to die off 
exponentially inside the medium, with a characteristic absorption length Labs- 

The known expression [2, 3, 12] of the absorption length can easily be recovered by 
means of the formalism exposed in the previous section, in the case of general anisotropic 
scattering and in the regime of weak absorption. We shall actually determine the extinction 
length of the more general problem defined by the coupled Q-dependent Schwarzschild- 
Milne equations (2.54). We look for slowly varying source functions of the form (2.58), 
along the lines of section 2.6. The main difference is that we have now wq = a ^ 1. It 
turns out that the estimates (2.60, 62) of the amplitudes ci^m still hold, whereas we obtain 

s^-Q^ + ^2_l^. (2.69) 

The Q-dependent extinction length in the presence of absorption reads therefore 

i^ext(Q) = - ~ ^ -YH (Q « 1, 1 - a « 1). (2.70) 

The usual absorption length is obtained by setting Q = in the above expression, namely 

/ pp* \ 1/2 

in agreement with refs. [2, 3, 12]. Another particular case is conservative scattering, in 
the absence of absorption, where we recover the result (2.63), namely 

i^ext(Q)~^~-- (2.72) 

This simple result holds for a general anisotropic scattering. It is a manifestation of 
the isotropic character of the long-distance diffusive behavior of the multiple scattering 
problem. We shall also come back to this point in section 5. 

3. LARGE INDEX MISMATCH REGIME 

In this section, we extend to the general case of anisotropic scattering the approach of 
ref. [15], which predicts the behavior of quantities in the regime where the optical indices 
no and ni of the diffusive medium and of the surroundings are very different from each 
other, i.e., when their ratio m = no/ni is either very small, or very large. As already 
pointed out in ref. [15], important simplifications occur in these regimes of a large index 
mismatch, where the Fresnel transmission coefficient of the boundaries of the medium is 
very small. To be more specific, radiation cannot enter the medium (respectively, leave 
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the medium) in the hmit m <^ 1 (respectively, m ^ 1), except at normal incidence. Ref. 
[12] already contains part of the results of this section. 

3.1. Diffuse reflection and transmission 

Along the lines of ref. [15], we evaluate the reflected and transmitted intensities in the 
large index mismatch regime by means of the following singular perturbative expansion of 
the Green's function Gs{t, ij,,t' , /j,'), defined by eq. (2.30). 

The starting point consists in noticing the identity 

ML{T,fX,T',fx')=R{-fx')MB{T,fX,-T',-fx') (3.1) 

between both kernels defined in eq. (2.27). Using -R(/x) — 1 — T{ij,), we can recast eq. 
(2.30) as 

Gs{r,ii,T',ii') = po{ii, ii')5{t - t') 

+ r^r" f ^^[MB(T-r",^,0,^")+MB(r + r",/i,0, V')] 

f3 2) 

- dr" ^T{^")Mb{t, a., -r", -^")Gs{t'\ r', /.'). 

In the limit of an infinitely strong index mismatch, i.e., for m = or m = oo, the 
transmission coefficient T(//) vanishes identically, so that the last integral of eq. (3.2), 
involving T(/x), is absent. It can be checked that the remaining terms only determine the 
Green's function up to an additive constant. This constant is only fixed by the action of 
the last integral, involving T(//), in eq. (3.2). It can therefore be expected to diverge as 
m — and m oo. 

In order to demonstrate this explicitly, we expand the Green's function according to 

Gs{t, /X, r', /x') = Cs + G'o(t, /X, r', /) + •••, (3-3) 

with the hypothesis that Cs diverges, whereas Go{t, fi,T' , /j,') remains finite, and the dots 
stand for terms which go to zero, as m ^ or m — > oo. The finite part Gq{t, n^r' , n') 
obeys the following equation 

Gq{t, II, t', ii') = poifi, ii')S(t - t') 

d/x 



+ J^^ dr" j ^^[Mb{t- t", II, 0, m") + Mb{t + t", yu, 0, -yu")] 



xGo{T",|^",T',|^') 



(3.4) 



-Cs I dr" I ^T{ui")MB{T,ui,-T",-ui"), 



-'0 
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together with the consistency condition 



dr dr' £^ ^ J^' '^T{^i')MB{T + r', 0, -//OColr', r", = 0, (3.5) 

derived along the hnes of ref. [15]. 

The constant Cs of the expansion (3.3) can be derived by integrating eq. (3.4) over 
the variables < r < oo and —1 < ji < 1. This yields 

Cs = ^, (3.6) 
where T is the mean flux transmission coefficient 



r = 2 ! iJ,T{n)dn 
Jo 



This quantity assumes its maximum T = 1 in the absence of any index mismatch, i.e., for 
m = 1, and it vanishes in both cases of a large index mismatch, according to 

(m< 1), 

(3.8) 

(m > 1). 



4m(m + 2) 



(m < 1), 



3(m + l)2 
4(2m + l) 



(3.7) 




The asymptotic behavior in the limits m ^ 1 and m ^ 1 of the quantities pertaining 
to reflection and transmission is immediately obtained by replacing in eqs. (2.33, 35) the 
Green's function Gs{t, ij,,t' , /j,') by its leading constant term Cs- We thus obtain 

4 4// 4//o/ib 

T"0 ~ Ti(/X) — , 7(/^a,/^b)~ ^ • (3.9) 



These predictions are identical to those derived in ref. [15], in the case of isotropic 
scattering. We have thus shown that the quantities which determine the diffuse reflected 
and transmitted intensities do not depend at all on the anisotropy of the cross-section in 
the large index mismatch limit. 

3.2. Enhanced backscattering cone 

The shape 7c (Q) of the cone of enhanced backscattering for a normal incidence can 
also be evaluated analytically in the regimes m ^ 1 or m ^ 1 of a large index mismatch. 
By inserting the estimates (3.9) into the general result (2.67), we flnd that the width of 
the cone is small in the large index mismatch regime, since it is proportional to the mean 
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transmission T. This observation suggests to consider the scahng regime where both Q 
and T are simultaneously small. 

In order to investigate this regime, we first recast the Q-dependent coupled Schwarz- 
schild-Milne equations (2.54) in a form similar to eq. (3.2). We then look for a solution 
of the form (2.58), where the inverse extinction length sq = Q is taken from eq. (2.63). 
We then proceed along the lines of section 3.1, integrating both sides of the coupled 
Schwarzschild-Milne equations over the variables < r < oo and —1 < fj, < 1. The inte- 
grals over r' which arc independent of the transmission T(/i) can be performed explicitly, 
whereas the Q-dcpcndcnce of the integrals involving T{p) can be neglected in the scaling 
regime. We thus obtain the following equations for the functions G^'^\Q, /j,) 

^G^-^HQ,/^) = Q5m,0 - Q £^yIo ^/^'^(/^')Pm(/^,/^')G("^^(Q,/^0 



-oo<fe<+oo 

with 



A(™''=)(g,/i) 



2 




Vl + Q2(l-/^2) V1 + + U {m>k). 

(3.11) 

The solution of eq. (3.10) in the scaling regime goes as follows. Along the lines of section 
2.7, we only keep the sectors m = and m = 1, and only the leading amplitude in the 
expansion (2.58) in each sector, namely we set G'^^\Q,fi) ~ co,o, G^^^^Q^fi) ~ Ci^i^/l — ix^. 
Inserting these expressions into eq. (3.10), all integrals can be performed, to leading order 
in Q. We thus obtain ci^i ^ — (Q/2)(r* — l)co,0) in agreement with eq. (2.60), while eq. 
(2.57) yields 7c(<5) ~ cq,q- After some algebra, we are left with the following scaling result 

ic{Q) ~ ^ ^Q^. (g«i,r«i). (3.12) 

T^~3~ 

The small-Q expansion of this prediction reads 

7c(g)«^^-^ + --- (3.13) 
The width of the cone therefore scales as 

3T 

Ag ^ — , (3.14) 

in agreement with the general formula (2.67), together with the results (3.9). 

The comment made at the end of section 3.1 still applies here. The scaling form of 
the enhanced backscattering cone in the regime of a large index mismatch does not depend 
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at all on the anisotropy of the scattering cross-section, apart from the simple power of r* 
already predicted by the diffusion approximation [31, 32]. 

4. VERY ANISOTROPIC SCATTERING 

In this section we investigate the regime where the scattering cross-section is very 
anisotropic, i.e., strongly peaked in a narrow cone of width ©rms 1 around the forward 
direction, with ©rms = (®^)- We thus have 1 — wi <^ 1, so that the anisotropy 

parameter r* ~ '^/^rms "^^^J large. The wavelength Aq of radiation in the medium, 
the scattering mean free path £, and the transport mean free path £* can therefore be 
considered as three independent length scales (Aq <^ £ <^ i*), besides other characteristic 
lengths, such as the sample thickness L, and possibly the absorption length Labs- 

The interest in this very anisotropic regime is twofold. First, we shall show that the 
radiative transfer problem in the absence of internal reflections is exactly solvable in this 
regime, just as it is for isotropic scattering, whereas the intermediate situation of a general 
anisotropy can only be treated numerically. Second, the dependence of physical quantities 
on anisotropy can be expected to yield the largest effects in the very anisotropic regime. 
We shall come back to this point in section 5. 

4.1. The example of large spheres 

We first recall how very anisotropic scattering can be realised experimentally. We 
consider the multiple scattering of light by large dielectric spheres, with radius a much 
larger than the wavelength Aq in the medium, i.e., with scale parameter koa ^ 1, so that 
geometrical optics can be used. Furthermore we assume that the optical index ns of the 
spheres is very close to the mean index uq of the medium, i.e., 

— = ms = l + (55 (|55|«1). (4.1) 
no 

The study of the scattering cross-section of electromagnetic waves by dielectric spheres 
is an old classical subject. The full solution was first derived by Mie in 1908. Ref. [33] 
provides an extensive overview of this field. The regime koa » 1 and \5s\ <^ 1 still has to be 
split into several subcases, according to the value of the combination \5s\koa. This can be 
understood as follows. In the framework of geometrical optics we can distinguish between 
the diffracted light, which is outgoing within an angle ©diffr ^ ^/{koa), independent of ds, 
and the refracted light, which is outgoing within an angle ©refr ~ independent of koa. 

From now on we concentrate our attention on the regime l^^l -C 1, koa ^ 1, and 
\5s\koa ^ 1. In this regime we have ©difir <S ©refr ^ 1- The cross-sections associated 
with each of the above processes asymptotically reads cTdiffr ~ Crefr ~ ^ra^, so that the 
total elastic cross-section a ~ 27ra^ is twice the geometrical one. This is the extinction 
paradox. We make the following approximations. We neglect the diffraction phenomenon 
by setting ©diffr = 0. We treat the refracted light according to the laws of geometrical 
optics, as illustrated in Figure 1. We neglect all the rays which are reflected at least once 
at the surface of the sphere, so that we only have to consider the refracted ray drawn on 
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the Figure. The corresponding phase function reads 



, . 4 sin /3 cos /3 

Prefr(0) = 



d/3 



de 



(4.2) 



sin O 

with the notations of Figure 1, which also imply 

e ^ -26stan/3. (4.3) 
Eqs. (4.2, 3) lead to the Lorentzian-squared scaling form [33] 

The cross-section is thus strongly peaked in the forward direction, as anticipated. 

We now turn to the evaluation of the coefficient wi^reti corresponding to refracted 
light. The integral (0^)refr = 6)^Prefr(6)) ©dO, with the phase function of eq. (4.4), 
is logarithmically divergent. A more careful treatment is therefore needed, which consists 
in using directly eq. (4.2), choosing u — sin^ jS as integration variable. We thus obtain 

1 



t^l.refr = ^ + ^ / udu^ {1 - w) (m| - u), (4.5) 



with Umax being the smaller of both numbers 1 and m|. The integral in eq. (4.5) can be 
performed explicitly for both ms < 1 and ms > 1. In both cases we obtain 

Wl,refr~ l-25|ln^, with A = 26-^/2 (l^sKl)- (4.6) 

The anisotropy parameter can be derived from the above estimate, not forgetting that 
CdifFr ~ Crefr ~ 0-/2 implies 1 — ■cc7i = (1 — zui^retr) f'^- To leading order for \5s\ <S 1, this 
yields 

£* 1 

= 7 - —nr- ("^ 

It is therefore evident that large spheres with a weak dielectric contrast provide a phys- 
ical instance of very anisotropic scattering, as described in the beginning of this section, 
with the refraction angle 0refr ^ -C 1 playing the role of 6rms) up to a logarithmic 
correction. We shall come back to this last point in section 5. 

4.2. Scaling limit of the Schwarzschild-Milne equation 

We now show that the general formalism of section 2 undergoes important simpli- 
fications in the regime of very anisotropic scattering. These will allow us to solve the 
Schwarzschild-Milne equation in the absence of index mismatch at the interface. This 
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analytical solution, to be described in sections 4.3 and 4.4, therefore shows that both lim- 
iting cases of isotropic scattering and of very anisotropic scattering are nearly on the same 
footing as far as the existence of exact results is concerned. 

The simplifications in the regime of very anisotropic scattering can be understood in 
physical terms as follows. Consider the random sequence of scattering events experienced 
by a light ray. At every scattering event, the direction n of the ray is only modified by a 
slight amount of order ©rms, with ©^ms ~ 2/t*. The extremity of the vector n therefore 
performs a Brownian motion on the unit sphere, with a small angular diffusivity e ~ ©rms- 
A similar picture has been used for long in the context of semifiexible polymer chains: the 
so-called Kratky-Porod theory models polymers as continuous persistent walks, whose unit 
tangent vectors obey a diffusion equation on the sphere (see ref . [34] for a review) . 

In more quantitative terms, the relationship (2.7) between the n-dependence of the 
specific intensity /(r, n) and the source function r(r, n) takes the form 

r(r,n) « (l + £A)/(r,n), (4.8) 

where A is the Legendre operator (Laplace operator on the unit sphere) 

d d'^ 1 

A = cot^— + ^ + -^^. (4.9) 

As a consequence, the integral operators T)^'^\ introduced in eq. (2.15), simplify to the 
following differential operators 

pM-i + eA^"*), (4.10) 

where 

A^^) = cot ^— + ^ - = 1 - _ 2 _ 4.11 

o9 89^ sm 9 o/j,^ d/i 1- 

is the Legendre operator in the sector defined by the azimuthal integer m. 

The angular diffusivity e is then determined by observing that Pi (//) = is an eigen- 
function of A*^'^-', with eigenvalue (—2), and of T)^'^\ with eigenvalue wi, by virtue of eq. 
(2.16). We thus identify wi = 1 — 2£, hence, using eq. (2.20), 

1 £ ©2 
2t* 2^* 4 ' ^ ^ 

in agreement with the above heuristic estimate. 

The radiative transfer equations (2.12) thus assume the differential form 

2T>-^/("^n^, A*) = A("^)/("^)(r, /x), (4.13) 
dr 

which demonstrates that quantities vary over length scales of order r ~ r*, i.e., z £*. In 
the next two sections we present our exact analytical predictions concerning observables 
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pertaining to optically thick slabs, in the regime of very anisotropic scattering, in the 
absence of internal reflections. The analysis roughly follows the presentation of section III 
of ref. [15], devoted to the case of isotropic scattering. 

4.3. Exact treatment in the absence of internal reflections: the homogeneous 
case 

This section is devoted to an analytic determination of the solution Th{t.,ij,) of the 
homogeneous Schwarzschild- Milne equation in the regime of very anisotropic scattering, 
and of the associated quantities of interest, tq and ti(//). 

It is convenient to consider the Laplace transform gH{s,fJ,) of Th{t,ij,), defined as 
follows 



I'OO 

9H{s,fJ,) = / dTTnir, -lJ,)e'^ (Res < 0). 
Jo 

The Schwarzschild- Milne equation (2.26, 28) can be recast as 



(4.14) 



1 + SjJ, 

9h{s,ii) -T*Ti{ll)/3 
1 + S/I 



(-l</x<0), 
(0</i<l), 



(4.15) 



since we have 5f//(— 1/^, /i) = t*ti(//)/3 for < < 1, as a consequence of eqs. (2.35, 39). 

We now expand eq. (4.15), using the expression (4.10, 12) of the operator V^^\ We 
introduce the rescaled Laplace variable 



a = 2sT*, 



as well as a new unknown function hn (cr, /x) , defined as follows 



(4.16) 



9h{s,ij,) = 



4r*2(l + sfj,)hH{(T, fJ.) (-K /X < 0), 

4T*2(l + S/x)/iH((T,/x)+T*Ti(/x)/3 (0 < < 1). 



(4.17) 



We assume that ti(//) vanishes faster than linearly in // for // ^ 0. This hypothesis will 
be checked a posteriori, as it will turn out that ti(//) ~ //^/^ for — > (see eq. (4.55)). 
The function hnio'ilj) is then continuously differentiable as a function of and it obeys 
the following equation 



(AW-a/.)/.^(a,/x) = | 



(-K/KO), 
Ti(//)/6 (0<//<l). 



(4.18) 



By means of the change of variable 



// = tanha;, (4-19) 
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eq. (4.18) can be recast in the form of the foUowing inhomogeneous Schrodinger equation 
on the real a;-hne 

- -^(-)''«(-. = { °v(x)K.)/6 > oj: <*-20) 
In this formula the accents denote differentiations with respect to x. The potential 

V{x) = tanhx(l - tanh^ x) (4.21) 
is an odd function of x, such that V{x)dx = /id/j,, and we have set 

Ti{ii) ^ liu{x) (0</i<l, x>0). (4.22) 

Equation (4.20) is a self-consistent equation for the two unknown functions huicr^x) 
and f^x). Analyticity properties in the complex cr- variable turn out to allow for an exact 
analytical solution of this equation. 

We introduce a basis of two elementary solutions {'Ui(cr, x), W2(<7, a;)} of the (homoge- 
neous) Schrodinger equation 

u"{a, x) - aV{x)u{a, x) = 0, (4.23) 

with the following asymptotic behavior 

ui{a,x)^l, U2{cr,x)^x (a; — > — oo), (4.24) 

up to exponentially small corrections. Similarly we introduce the basis {vi (cr, x), V2{(t, x)}, 
with the asymptotic behavior 

vi(cr, a;)~l, V2{o;x)^x (a; — > oo). (4.25) 

The Schrodinger equation (4.23) admits solutions with the above boundary conditions, 
since the potential V{x) vanishes exponentially as x ^ ±oo. The four above solutions are 
entire functions of a, i.e., they are analytic in the whole cr-plane. Furthermore they are 
related by the following identities 

vi{a,x) = ui{-a,-x), V2{(t, x) = -U2{-(t, -x), (4.26) 

because the potential V{x) is odd. 

We recall that, if u{x) and v{x) are any two solutions of the Schrodinger equation 
(4.23), with the same cr, their Wronskian 

W{u, v} = u{x)v'{x) - u'{x)v{x) (4.27) 
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is independent of x. The boundary conditions (4.24, 25) imply that both bases of functions 
have unit Wronskian, namely 



W{ui,U2} = W{vi,V2} = l. (4.28) 

For generic values of the parameter cr, both bases of solutions are related by a 2 x 2 
transfer matrix of the form 

V2j- \H{a) F{-a)J \u2 J ' ^^"^^^ 

The three functions which enter eq. (4.29) are entire functions of cr; the determinant of 
the matrix is F{a)F{—a) — G{a)H{a) = 1; as a consequence of eq. (4.26), G{a) and H{a) 
are even functions of cr. 

The above functions also govern the non-trivial asymptotic behavior of the bases of 
solutions 

ui{a,x)^ F{-a)-G{a)x, U2{a,x) ^ -H{a) + F{a)x {x ^ oo), 
vi{a,x) !v F{a) + G{a)x, V2{o;x) !v H{a) + F{—a)x {x^—oo), 

as well as their mixed Wronskians 

W{ui,v,} = G{a), W{ui,V2}=F{-a), W{u2,v^} = -F{a), W{u2,V2} = -H{a). 

(4.31) 

The first terms of the Taylor expansion of wi(cr, x) and W2(cr, x) around cr = read 

ui(a, x) = 1 — —{1 + tanh x) + ■ ■ ■ , 

^ rx X (4-32) 

tt2(cr, x) = X + a — tanhx) + ln(2 coshx) j + • • • 

We have also determined the terms of order cr^, which are too lengthy expressions to be 
reported here. They imply 

F(ct) = 1+(J+(jV2 + ---, ^((j) = ctV3+---, if((j) = (7/6-7rV36)(j2 + --- (4.33) 

On the other hand, large values of the complex parameter cr correspond to the semi- 
classical regime for the Schrodinger equation (4.23), where the behavior of the functions 
ui{a,x) and U2(cr, x) can be derived by means of a W.K.B.-like approximation. This 
regime is analyzed in detail in Appendix B. Let us mention that the wavefunctions display 
oscillations when either cr > and a; < or cr < and x > 0, whereas they are growing or 
decaying exponentially in the other two cases. 

The function G{a) deserves some more attention, since it will play a central role in the 
following. G{a) can be viewed as the functional determinant of the Schrodinger equation 
(4.23), in the following sense. Assume cr is such that G{a) = 0. We have then 

vi{a,x) = F{a)ui{a,x), ui{a,x) = F{-a)vi{a,x), F{a)F{-a) = 1. (4.34) 
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In other words, for such a cr, the Schrodinger equation (4.23) has a bounded solution 
over the whole real line. Throughout the following, using slightly improper terms, such 
a value of a is called an eigenvalue of the Schrodinger equation, and the corresponding 
function fi(cr, x) is referred to as the associated eigenf unction. The semi- classical analysis 
of Appendix B demonstrates that there is an infinite sequence of real eigenvalues. We label 
them by an integer — oo < n < oo, so that Un > for n > 1, uo = 0, and a-n = —an- The 
associated eigenfunctions vi{an,x) are orthogonal with respect to the following indefinite 
metric 

/"OO 

vi{am,x)vi{an,x)V{x)dx = NnSm,n (-00 < m, n < oo) . (4.35) 



The squared norms obey the symmetry property 



(4.36) 



as a consequence of eq. (4.34). The eigenfunction Vi{Q,x) = 1 associated with the zero 
mode ctq = is peculiar, since its squared norm reads Nq — Q. As a consequence, the basis 
of eigenfunctions {vi(cr„, x), n ^ 0} only spans the set of bounded functions f{x) on the 
real x-line such that 

/■oo 

f{x)V{x)dx = 0. (4.37) 



For such functions, we have 



with 



/(^) = XI ^nVlicTn, X), (4.38) 



1 

An7_c 



vi{an,x)f{x)V{x)dx. (4.39) 

-oo 

The determination of the contribution of the zero mode to functions f{x) which do not 
obey the condition (4.37), such as e.g. a constant, requires more care. An elegant way of 
dealing with this problem consists in introducing a deformation parameter k,, as we shall 
see in section 4.4. 

It is advantageous to factor the entire function G{a) as follows 

Gia) = ^P{a)P{-a), (4.40) 

with the notation 

pi<r) = n f 1 + ^) • (4.41) 



n>l 



n 



The explicit solution of the inhomogcncous Schrodinger equation (4.20) now goes as 
follows. Since hu^cr^x) is a regular function of x in the a; — > — oo limit, it is proportional 
to wi(a", x) for a; < 0, namely 

hH{(T,x) = aH{o-)ui{a,x) (a; < 0). (4.42) 
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For X > we solve eq. (4.20) by varying the constants, namely we look for a solution of 
the form 

hnio'jx) = bH{o',x)vi{a,x) + CH{cr,x)v2{(J,x) (a; > 0). (4.43) 
The unknown constants hH{<y,x) and c//(cr, x) obey the requirements 

h'jj{a,x)vi{a,x) -\- c'jj{a,x)v2{(T-,x) = 0, (4.44a) 
b'jj{(T,x)v[{(7,x) + c'jj{(7,x)v2{<7.,x) = V{x)v{x)/Q, (4.44b) 

where the accents again denote differentiations with respect to x. Eq. (4.44a) is a con- 
straint imposed a priori, in order to break the large redundance of the representation 
(4.43); eq. (4.44b) is then a consequence of eq. (4.20). Eq. (4.44) can be solved for h'jj 
and c'jj, using eq. (4.28). Integrating the expressions thus obtained, we get 



poo 

bnicr, x) = 6//(cr, oo) + / V2{(t, y)iy{y)V{y)dy/6, 

J X 

/"OO 

CH{cr,x)^- vi{a,y)iy{y)V{y)dy/Q. 

J X 



(4.45) 



Indeed there cannot be a non-zero constant ch{o', oo), because this would correspond to 
an unacceptable singular solution of the form /iij(cr, a;) ~ a; ~ ln(l — /x) for //—>!. 

Both expressions (4.42) and (4.43) have to match at a; = 0, together with their first 
derivatives. These two conditions determine 6h(c, 0) and ch(c, 0), and some algebra then 
leads to the identity 

/"OO 

G{a)aH{cr) ^ -CH{cr,0) ^ vi{a,x)iy{x)V{x)dx/6. (4.46) 

Jo 

We can argue on eq. (4.46) as follows. Since qh^s^ij) is analytic in the half-plane 
Res < 0, auicf) has the same property for Reu < 0, so that the zeroes — cr„ of G{a) 
cannot be poles of anicr). They are therefore zeroes of Cif(cr, 0). On the other hand, the 
a^s are not zeroes of c//(o", 0), since the integral expression in eq. (4.46) is positive for 
cr > 0. Hence they are poles of anio')- The final step concerns the small-cr behavior of the 
quantities involved in eq. (4.46). The asymptotic behavior (2.29) of rij(T, /i) yields the 
following double-pole structure for its Laplace transform near s = 

gH{s,fi) = \- ^*^^'^''K oil) (s^O), (4.47) 

s 



which implies the following small-cr behavior of aH{<y) 

1 1 



aH(a) = ^ + ^-^ + 0(l) (cT^O). (4.48) 



We thus obtain finally 
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Eq. (4.49) can be considered as an explicit result. Indeed -P(cr), defined in eq. (4.41), is for 

all purposes a known function, since the eigenvalues can be determined numerically, es- 
sentially with arbitrary accuracy, via the partial- wave expansion procedure of Appendix A. 
Furthermore the semi-classical analysis of Appendix B determines the asymptotic behavior 
of the eigenvalues in the regime of large quantum numbers (n 3> 1). 

As a first consequence of the above results, we can determine the reduced thickness tq 
of a skin layer, i.e., the reduced extrapolation length, by comparing the small-u behavior 
of the expression (4.49) for anicr) with the expansion (4.48). We thus obtain 

To = 1 - 2 y — = 0.71821164. (4.50) 

n>l 

The series converges, since cr„ grows as n^, according to the semi-classical estimate (B.21). 
The number given in eq. (4.50), as well as all the subsequent ones, has been obtained by 
means of the partial- wave expansion described in Appendix A. This number gives an idea 
of the accuracy of this approach. The most significant numerical results are listed in Table 
2, together with their counterparts in the case of isotropic scattering. 

Second, the determination of Ti{p) goes as follows. We recall that this quantity is 
related to i^{x) by eq. (4.22). Eqs. (4.46, 49) yield 

/■oo 

/ vi{a,x)iy{x)V{x)dx = 2P{a). (4.51) 
Jo 



The small-cr expansion (4.32) of ui{a, x), together with eq. (4.26), allows us to recover 
the sum rules (2.38, 42) in the following form 

/•oo nOO 

/ i^{x)V{x)dx = 2, / i^{x)tanh.xV{x)dx = 2To. (4.52) 
Jo Jo 



On the other hand, the semi- classical estimates (B.18, 20) of vi{a,x) and -P(cr) for 
large values of cr = K'^ imply 

POO 

/ u{x)Ai{K^/^x)xdx^{6/7ry/^K-^/^ (K>1). (4.53) 
Jo 

This estimate yields, by means of eq. (B.25) for s = 3/2, the small-x behavior of i'{x), i.e., 

i^{x) ^ 6{2x/7rf/^ {x < 1). (4.54) 



We thus obtain the following universal scaling behavior of the Ti-function in the case 
of very anisotropic scattering 

Ti(//)Ri6(2/7r)^/V/' (4.55) 
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This novel result is in contrast with the linear behavior ri(yLt) ~ /i-\/3 observed for isotropic 
scattering. The law (4.55) also confirms the hypothesis made in the beginning of this 
section, namely that ri(//) vanishes faster than linearly in fi. 

Finally, we can extract the full functions i'{x) and ti(//) from eq. (4.51) by means of 
the inversion formula (4.38, 39). We obtain 

iy{x) = 3(to + tanhx) + 2 ^ ^^^vi{an, x), (4.56) 



n>l ' 



I.e., 



lM = 3^ro + f^) + 2j2 ^P^v,{an, argtanh/.). (4.57) 

^ n>l 

The semi-classical analysis of Appendix B implies that the contribution of the n- 
th mode to the results (4.56, 57) falls off exponentially with n for x > 0, according to 
exp [ — Kn{l — /(x))]. This exponential convergence disappears as a: — 0, where eqs. 
(4.54, 55) apply. The explicit terms in front of the sums in eqs. (4.56, 57), corresponding 
to the contribution of the zero mode (n = 0), have been anticipated from results to be 
derived in section 4.4. 

Figure 2 shows a plot of the function Ti(/i), for both isotropic scattering [15] and very 
anisotropic scattering (eq. (4.57)). The maximal values ti(1) for both cases are given in 
Table 2. The difference between both limiting cases is remarkably small. 

4.4. Exact treatment in the absence of internal reflections: the inhomogeneous 
case 

In this section we derive analytical expressions in the regime of very anisotropic 
scattering of the special solution Ts{t, /i, /la), with its by-product the bistatic coefficient 

In analogy with eq. (4.14), we define the Laplace transform of the source function as 
follows 

I'OO 

gs{s, II, iia) = dTTs{T,-fi,fia)e'^ (Res<0). (4.58) 

Jo 

The Schwarzschild- Milne equation (2.26, 28) can be recast as 



9s{s,II,Ha) 



^aPo(-/^,/^a) ^ p(0) 
1 - SHa 



1 + S/i 
gs{s,II,lia) - 'fill, l^a)' 



(-l<//<0), 

{0<fi<l), 



(4.59) 



I + S/J, 

since we have gs{—'^/n-,lJ-,lJ-a) = lil^a) for < < 1, as a consequence of eq. (2.33). 

In analogy with section 4.3, we expand eq. (4.59), using eqs. (4.10, 12). We use the 
rescaled variable cr of eq. (4.16), and we introduce a new unknown function hsicr, fi, fia), 
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defined as follows 



^ r« /, /, \ _ / 2t*(1 + s^)hsi(7, //, (-K /i < 0), . . 

gs[s, - I 2r*(l + + 7(a*, A^a) (0 < < 1). ^^"^^^ 

We assume that 7(/ia5A*6) vanishes faster than linearly as //o or — > 0. Using the change 
of variable (4.19), we are again left with an inhomogeneous Schrodinger equation, namely 

'•sf''- -.-')={;!n*wr.j (^>o): 

where we have set 

7(//,//a) = iiiiaP{x,Xa) (yU = tauhx > 0, = tauhxo > 0). (4.62) 

Now, in order to deal with the problem of the zero modes of the Schrodinger equation 
(4.23), we introduce a continuous deformation parameter k > as follows. We consider 
the deformed Schrodinger equation 

u"{a, x) — {aV{x) + K'^W{x))u{a,x) = 0, (4.63) 

with 

W(x) = (1 - tanh^ xf. (4.64) 

The eigenvalues with label n ^ acquire a regular ^-dependence of the form 

anin) = -a^n{n) ^(Tn + 0{K^), (4.65) 

as well as the associated eigenfunctions vi{k, a"„(re),a;), whereas the double degeneracy of 
the zero mode ctq = is lifted into the following two exact eigenvalues and eigenfunctions 
ofeq. (4.63) 

a±o{K) = ±2k, vi{k, o^±o{K), x) = exp ( ± k{1 — tanhx)) = exp ( ± k{1 — //)), (4.66) 
with squared norms 

/ X , e"^^"^ /sinh2«: , ^ \ , . 

N±o{k) = ± — cosh 2k . (4.67) 

The introduction of labels ±0 is consistent with setting n = in formulas such as eq. 
(4.36) or (4.65), rather than with the standard arithmetics of integers! 

For any non-zero the set of eigenfunctions {fi(K, cr^(«:), a;)}, where the label n runs 
over the non-zero algebraic integers (n ^ 0) plus both values n — ±0, now spans the whole 
space of bounded functions j{x) on the real line. The difficulty of the constraint (4.37), 
due to the vanishing norm of the zero mode at k = 0, is thus cured in a natural way. 
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The mixed Wronskian G{n, a) = W{ui,vi} can still be factorised over its zeroes, in 
analogy with eq. (4.40), namely 

G{k, a) = G{k, 0)R{k, a)R{K, -a), (4.68) 

with 

= I ('"^) 

The prefactor G{k, 0) of eq. (4.68) is a non-trivial function of k. Indeed this quantity 
can be viewed as the functional determinant of the Schrodinger equation (4.63) with cr = 0, 
namely 

u"{a, x) - K'^W{x)u{a, x) = 0. (4.70) 

Eq. (4.70) is equivalent to the spheroidal equation studied by Meixner and Schafke [35]. 
Some properties of this equation have been studied in detail [36], in an investigation of 
nematic phases of semiflexible polymer chains. The occurrence of eq. (4.70) in that context 
is related to the Kratky-Porod description of persistent chains, mentioned in section 4.2. 

Eq. (4.70) has a discrete spectrum of imaginary eigenvalues of the form k = ±i^„ 
(n > 1), as shown by the semi-classical analysis of Appendix B. On the other hand, the 
regularity of G(k, cr) at k = implies the small-K behavior G(k, 0) 4k;^/3, hence 

^(«'0) = ^n(l + ^)- (4-71) 

„>i V ?n/ 



(4.69) 



The solution of the K-dependent deformed inhomogeneous Schrodinger equation (4.61) 
then follows the lines of section 4.3. The particular values x — —Xa < and x — define 
three sectors, in which we look for a solution of the following form 

{as{K, cr, Xa)ui{K, cr, x) [X < -Xa), 

bsiK, cr, Xa)ui{K, cr, x) + Cs{k, a, Xa)u2{K, Cr, x) {-Xa < X <0), 

ds{K, cr, X, Xa)vi{K, cr, x) + es{n, CF, X, Xa)v2{K,, cr, a;) {x > 0). 

(4.72) 

The constants which enter the last of these expressions obey the conditions 

d'g{K, cr, X, Xa)vi{K, cr, x) + e'g{K, cr, x, Xa)v2{K, cr, x) = 0, 

d's{K,a,X,Xa)v[{K,a,x) + e's{n, a,X,Xa)v'2{n, a,x) = IJ,aV{x)p{K,X,Xa), 

hence 

/•oo 

ds{K,a,x,Xa) = ds{K,a,oo,Xa) + Ha / V2{K,a,y)p{K,y,Xa)V{y)dy, 

J X 

poo 

es{K, cr, X, Xa ) = -Pa / vi{k, cr, y)p{K, y, Xa)V{y)dy. 

J X 
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(4.73) 



(4.74) 



On the other hand, the matching of the solution (4.72) at a; = —Xa yields 

Cs{K,a,Xa) = -2iiaUi{n,a,Xa), 

whereas its matching at a; = leads to 

ds{K, cr, 0, Xa) = F{k, -a)bs{K, a, Xa) ~ H{k, a)cs{K, cr, Xa), 
-es{n, a, 0, Xa) = G{k, cr)bs{n, a, Xa) - F{k, cr)cs{n, a, Xa), 



(4.75) 



(4.76) 



with the notations (4.29). We are thus left with 

G{k, a)as{n, cr, Xa) = -2ij,aVi{K, a, -Xa) - es{n, cr, 0, Xa). (4.77) 

We can now follow the approach used on eq. (4.46). Since a5(«:, cr, Xa) is holomorphic 
in the half-plane Re a < 0, the zeroes — a"^(«:) for n > of G(k, cr) cannot be poles of 
as{n, cr, Xa)- Hence they are zeroes of the right-hand side of eq. (4.77). 

We are therefore left with the problem of finding es{K, o", 0, Xa), an entire function of 
a, from the knowledge of its values on the sequence of points cr = —crn{><') {n > 0), together 
with a natural assumption of minimal growth at infinity compatible with these data. This 
is a generalisation to an entire function of the problem of finding a polynomial Q{z) with 
minimal degree, knowing its values at N points, namely Q{zn) = Qn for 1 < n < N. It is 
useful to view the z^s as the zeroes of the normalised polynomial 

Piz)= n i^-^n). (4.78) 

l<n<N 

The solution Q{z) with minimal degree (generically N — 1) is given by the following La- 
grange interpolation formula 

l<n<N l<m^n<N l<n<JV ^^^^ ^^^^ 

Extending eq. (4.79) to the present case of an infinite sequence of data for an unknown 
entire function, we obtain 

-es{n, a, 0, Xa) = 2fiaR{n, (t) V . . ^/''^w 10"/^ u — T^' '^^•^°) 

^ (C7 + C7„(K))(di?/dc7)(K, -C7„(k)) 

or equivalently, using a generalisation of the identity (C.6) to k 7^ 0, together with the 
definition (4.68), 

/ n ^ on/ R{K,an{K))vi{K,an{K),Xa) , . 

-es{K, a, 0, Xa) = -2^aG{K, 0)R{k, a) > — ■ . (4.81) 

^ [a + an[K))Nn[K) 
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Finally, we can derive an explicit expression for p{K,x,Xa), by means of an inversion 
formula analogous to eqs. (4.38, 39), namely 

/ N o^/ n\ ^ R{K,am{n))R{K,an{n)) Vi{K,am{n),x) Vi{K,an{K,),Xa) 
p{K,X,Xa) = -2G{k,0) > ^ —— 

m,n>0 



E 



vi{k, an{n),x)vi{K, an{n), -Xa) + vi{k, an{n),Xa)vi{K, an{n), -x) 



(4.82) 

This central result is manifestly symmetric under the exchange of x and Xa, as it should. 
This symmetry property is a valuable check of the whole approach, since both arguments 
X and Xa have played uneven roles throughout the derivation. 

We are now able to take the physical k — > limit of the above results. In this regime 
eq. (4.81) can be recast as 



\ n>l ^ 



-esia, 0, Xa) = 2piaP{cT) I 1 - ^ > : I • (4-83) 



First, we are now able to complete the proof of the anticipated result (4.56, 57), by inserting 
eq. (4.83) into eq. (4.77), expanding the latter equation for k — > to the first non-trivial 
order as u — > 0, and comparing the result with the estimate 

as{a,xa) ^ -^-^ (a^O). (4.84) 
cr 

Second, the small-cr expansion of eq. (4.83) allows us to recover the sum rules (2.37, 41) 
in the following form 

poo poo 

/ p{x,Xa)V{x)dx = 2, / p(a;, a^a) tanha;y(a;)da; = 2z>'(a;a)/3 — 2 tanhxa. (4.85) 
Jo Jo 

On the other hand, for large values of cr = K^, we can use the semi-classical estimates 
(B.18, 20, 22), setting cr„ = /c^, in order to transform eqs. (4.74, 83) into 

/•OO poo u2/3 

/ p{x,Xa)Ai{K^/^x)xdx^{2/n)K^/^ Ai{k'^/^Xa)dk. (4.86) 

Jo Jo «^ + -f^ 

This estimate shows that xxap{x, Xa) is a homogeneous function of its arguments with 
degree zero, when both of them are small, i.e., xxap{x, Xa) — g{x/xa)- The rescaling of 
eq. (4.86) according to 2; = K'^/^x = k'^^^Xa then yields, by means of a mere identification 
of both integrands, using eq. (B.25), the expression g(z) = (3/77)2:^/^/(2;^ + 1), implying 
the following scaling behavior 



7r{xl + xf) 



p{Xa, Xt) ~ _C 3X i^a, Xb « 1), (4.87) 
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or equivalently 

7(//a, l^b) ~ ^[^ffl il^a, l^b « 1). (4.88) 

7r(/io + fJ'b) 

This novel result is in contrast with the rational behavior 7(/Ua5 A^fc) ~ fJ-afJ-b/iPa + A^fo) 
in the case of isotropic scattering. The law (4.88) confirms the hypothesis made in the 
beginning of this section, namely that 7(/Xa, /x^) vanishes faster than linearly in either of its 
arguments. It is also worth noticing that the scaling form (4.88) of the bistatic coefficient 
saturates the sum rule (2.37). 

The full expression of the bistatic coefficient 7(//a) A*6) is obtained by taking the k — > 
limit of eq. (4.82), using the definition (4.62). The modes m, n = yield divergent 
contributions as k 0, which cancel out as they should, as well as finite parts, so that we 
are left with the result 



lii^a.^lb) 



l^al^b 2 n>l 



^ ^[viicTn, Xa)vi{an, -Xb) + Vi{an, -Xa)vi{an, Xb)] 



n>l 



(4.89) 



2 v-^ O'mO'n -P(o'm) -P(o'n) / x / n 

o TT^rZ AT ir^M(^m,Xa)Vi{an,Xb), 



with Ha = tanhxa, lib = tanhx^. 

The maximal value of the bistatic coefficient, which yields an absolute prediction for 
the diffuse reflected intensity in the normal direction, reads 7(1, 1) = 4.889703. This 
number is some 15% above the corresponding one in the case of isotropic scattering (see 
Table 2). 

4.5. Extinction lengths of azimuthal excitations 

Up to this point, we have mostly investigated quantities with cylindrical symmetry 
around the normal to the slab, pertaining thus to the m = sector of the azimuthal 
decomposition (2.8). 

We now want to consider briefly the other values of the azimuthal integer m, in the 
regime of very anisotropic scattering. As already mentioned in section 2, all the sectors 
contribute e.g. to the reflected intensity, except if the incident beam is normal to the 
sample, or, more generally, has itself cylindrical symmetry. The situation is different in 
transmission through thick slabs, to which only the sector m = contributes. The reason 
for this is that the intensity in the other sectors is exponentially damped inside the sample, 
namely 

/(-)~exp(-z/Lext^"^)), (4.90) 
where Lext*-"^^ is the extinction length of the azimuthal excitations in the sector m. 

It is the purpose of this section to determine these lengths in the limit of very 
anisotropic scattering. The Legendre operator in the sector deflned by the azimuthal 
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integer m is given by eq. (4.11). As a consequence, and along the lines of sections 4.3 and 
4.4, we are led to study the Schrodinger equation 

u"{a, x) - (m^ + aV{x))u{a, x) = 0. (4.91) 

The corresponding extinction length is given by 

9/* 

= -^y (4.92) 

where the smallest positive eigenvalue of eq. (4.91). This is indeed the location 

of the first singularity of the Laplace transform of the intensity. 

For large values of the azimuthal integer m, we can use the semi-classical analysis 
developed in Appendix B. The Sommerfeld quantisation formula associated with eq. (4.91) 
reads 

/ {-aV{x)-m^y^^dx= / (a-^- ^ ] d/x ^ (n + l/2)7r. (4.93) 

The smallest eigenvalue (Jo^"^-' corresponds to setting n = in the above formula. In first 
approximation we express that the argument of the square-root inside the /x-integral in eq. 
(4.93) has zero as its maximal value, so that = We thus obtain uo^"^^ ^ 3-\/3m^/2. 
In second approximation we expand the integrand around its maximum, which takes place 
for fj, a; \/3/3. We obtain after some algebra the following next-to- leading order estimate 

ao^"^^ ^ (3\/3/2)(m2 + mV2), (4.94) 



hence 

4/* 

LeJ"^^ ^ ^ (m > 1). (4.95) 

This semi-classical estimate gives accurate numbers of the whole spectrum of ex- 
tinction lengths, down to the largest one, Lext*^^"*- Indeed eq. (4.95) yields Lext*^''^'* ~ 
0.318861 £*, whereas the exact numerical value reads Lext^^^ = 0.2829169 £*. 

For a large but finite anisotropy, the lengths Lext*^'"-* follow the universal law (4.95) 
only for m < m* ~ -\/t* ~ l/6rms- For larger values of the azimuthal number, the Lgxt*^"^-* 
become non-universal numbers of order £. This crossover is expected on physical grounds. 
Indeed large azimuthal numbers m ^ m* correspond to an angular resolution SQ <^ ©rms) 
so that the details of the cross-section matter in this regime. 

5. DISCUSSION 

In this paper we have considered several aspects of multiple anisotropic scattering of 
scalar waves. We have considered the geometry of an optically thick slab, of thickness 
L ^ £*. Our main goal has been to investigate in a quantitative way the effects of the 
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anisotropy of the scattering cross-section, and of the internal reflections at the boundaries 
of the sample, due to an optical index mismatch. 

The general results derived in section 2 show that, in first approximation, quantities 
only depend on the anisotropy through the transport mean free path £*. This is especially 
the case for the angle-resolved transmission through a thick slab (2.47), for the thickness 
of a skin layer, zq — tqI* , and for the width (2.68) of the enhanced backscattering cone. 
The present work thus confirms on a firm basis that the scaling behavior of these quanti- 
ties is qualitatively explained within the diffusion approximation, which amounts to only 
considering the long-distance diffusive character of the propagation of radiation in a turbid 
medium. The scaling law in l/£* of the width of the cone has been derived long ago within 
the diffusion approximation [31, 32]. The results of section 2.7 concerning the dependence 
of the extinction length with respect to Q and a can be directly compared with the pre- 
diction of the diffusion approximation [12]. Within this framework, all extinction effects, 
provided they are small enough, can be coded in a single parameter, namely the mass M. 
such that 

= q^ + iVL+^, (5.1) 

abs 

where q is the transverse wavevector. Labs is the absorption length, and Vt — (w — c<;')/-Dphys 
represents the properly dimensioned contribution of a small frequency shift between the 
advanced and the retarded amplitude propagators which build up the diffuson. The inverse 
extinction length is then equal to the real part of the complex mass Ai. Our results fully 
agree with eq. (5.1), with A4 ~ so/i, q — Q/^, and Labs as in eq. (2.71). 

We have then investigated in detail to what extent observable quantities are univer- 
sally described by their explicit dependence on £* recalled above, and to what extent they 
still depend on details of the scattering cross-section mechanism. As recalled in the Intro- 
duction, this question is beyond the scope of the diffusion approximation, and requires a 
careful treatment using the radiative transfer theory, at least in the regime £ 3> Aq. For the 
diffuse reflected or transmitted intensity, and for the width of the enhanced backscattering 
cone, the detailed structure of the scattering mechanism only contributes a small effect, 
entirely contained in prefactors of the laws mentioned above, such as the constant tq, or 
the functions ti(//) and 

Two regimes of interest allow for more quantitative results: 

(i) The regime of a large index mismatch, where the boundaries of the sample almost act 
as perfect mirrors, is considered in section 3. Our results (3.9, 12) are identical to those 
derived in ref. [15], in the case of isotropic scattering. Therefore the quantities we have 
considered do not depend at all on the scattering cross-section in this regime. This can 
be understood as follows. Since the thickness Zq ~ 4£*/(3T) of a skin layer is very large, 
the radiation undergoes many scattering events near the boundaries before it leaves the 
medium, so that the details of every single scattering event are washed out. 

(ii) In the absence of internal reflections, we have considered in detail the regime of very 
anisotropic scattering. In section 4 we have presented an exact analytical treatment of the 
radiative transfer problem in this regime. We have obtained the results (4.50, 57, 89) which 
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determine the diffuse reflected and transmitted light for a thick slab. These results are given 
in terms of the eigenvalues and eigenfunctions of onc-dimcnsional Schrodinger equations, 
which are accessible both numerically, via the partial-wave expansion of Appendix A, 
and analytically in the limit of large quantum numbers, via the semi-classical analysis of 
Appendix B. It is a priori possible to extend this exact treatment to the scaling behavior 
of the shape of the enhanced backscattering cone. The general structure of the equations 
to be solved shows that we have 7c(Q) ~ F{Qt*), in a whole scaling region defined by 
Q <^ 1 and r* ^ 1. By inserting the numerical values of Table 2 into the expansion 
(2.66), we get at once F(0) = 7(1, 1) = 4.889703 and F'(0) = -ti(1)73 = -8.80166. The 
exact determination of the full scaling function F would amount to solving a self-consistent 
inhomogcncous equation of the type (4.61), albeit with the full Legendre operator instead 
of one-dimensional second-order derivative. The wings of the cone, starting around values 
of Q of order unity, will depend on the details of the scattering cross-section, even in the 
regime of very anisotropic scattering. 

Our exact treatment of the radiative transfer problem in the very anisotropic regime, 
based on the expansion (4.8), is expected to be valid in the regime ©rms <S 1 of a broad 
universality class of phase functions. Although this universality class cannot be easily 
characterised, we can assert that it contains at least the phase functions scaling as 

p(e)^$(e/e™s), (5.2) 

such that the scaling function $ has a finite second moment. This restrictive definition 
docs not encompass a priori the Lorentzian-squared phase function (4.4), which has a 
logarithmically divergent second moment, as already mentioned in section 4.1. The same 
remark holds for the so-called Henyey-Greenstein phase function 

^^^^ = (1-2^ cos 6 + ^2)3/2' ^^-3) 

often used in numerical investigations [3, 17], for which the second-moment integral is 
linearly divergent. 

The discussion of the dependence of quantities on the details of the scattering mech- 
anism is summarised in Table 2, where we compare the numerical values of the dimen- 
sionless absolute prefactors of five characteristic quantities, for isotropic scattering and for 
very anisotropic scattering. The relative differences, shown in the last row, are very small 
in most cases. Some other quantities, such as the shape of the enhanced backscattering 
cone, or the spectrum of extinction lengths of the azimuthal excitations, exhibit universal 
behavior in the very anisotropic regime only in a limited range, corresponding to a low 
enough angular resolution (50 ^ ©rms); so that the details of the scattering cross-section 
do not matter. 

Finally, we can compare our universal results in the very anisotropic scattering regime, 
for some of the quantities listed in Table 2, with the outcomes of numerical approaches. Van 
de Hulst [3, 17] has investigated in a systematic way the dependence of various quantities 
on anisotropy, for several commonly used phenomenological forms of the phase function. 
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including especially the Henyey-Greenstein phase function (5.3). The data on the skin- 
layer thickness reported in ref. [3] show that, as a function of anisotropy, Tq varies from 
0.7104 (isotropic scattering) to 0.7150 (moderate anisotropy), passing a minimum of 0.7092 
(weak anisotropy). The trend shown by these data suggests that our universal value 
0.718211 is actually an absolute upper bound for tq. Numerical data concerning ri(l) 
is also available. Van de Hulst [17] has extrapolated two series of data, concerning the 
Henyey-Greenstein phase function (5.3), which admit a common limit for very anisotropic 
scattering {g ^ 1). According to the analysis of section 4, this limit reads in our language 
Ti(l)/4 = 1.284645, whereas ref. [17] gives the two slightly different estimates 1.273±0.002 
and 1.274 ± 0.007. The agreement is satisfactory, although it cannot be entirely excluded 
that the observed 0.8% relative difference can be a small but genuine nonuniversality effect. 
Indeed, as mentioned above, the Henyey-Greenstein phase function (5.3) might not belong 
to the universality class where our approach holds true. The same remark applies to a 
less complete set of data [17] concerning the intensity 7(1, 1) of reflected light at normal 
incidence. 
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Appendix A. Partial-wave expansions 

In this Appendix we describe a numerical algorithm based on a partial- wave expansion, 
that we have used to determine the eigenvalues and eigenfunctions of the Schrodinger 
equations involved in section 4. 

We first consider the Schrodinger equation (4.23). Going back to the //-variable, this 
equation reads 

(A(o)-a//K(a,//)=0, (A.l) 

where A*^*^) is the Legendre operator in the (/^-independent sector, defined in eq. (4.11). It 
is natural to expand the function vi {a, jj) in the Legendre polynomials 

^;i(c7,//) = J]a^(a)P^(//). (A.2) 

Indeed these polynomials are eigenfunctions of namely 

A(o)p^ = -£(£ + l)P^, (A.3) 
and the product iJ,Pe{iJ,) has the following expression 

(2£ + 1)/.P,(//) = {e+l)P^+^{|^)+eP^_^{|l), (A.4) 
so that eq. (A.l) amounts to the following three-term recursion relation 

£{£ + l)ai + a A^a,+i + = q. (A.5) 

When a is one of the eigenvalues eq. (A.5) has an acceptable solution {aeicr)}, 
decaying to zero for large £. The quantities needed in section 4 can then be evaluated as 
follows. The normalisation condition (4.25) becomes 

^a^((7„) = l, (A.6) 

since Pi{l) = 1. The squared norms of the eigenfunctions read 

= 4|: ^,^/,^^,'^^3) a,(cx>mK), (A.7) 

as a consequence of the normalisation of the Legendre polynomials 
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' ^iPMPM = (A.8) 
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Finally, the non-trivial mixed Wronskians F(±cr^) read 

F{-an) = = ^iK,/^ = -1) = (A.9) 

since = (-1)^. 

We now consider the K-dependent Schrodinger equation (4.64), namely 

{A^°^ -a/i + K^ili^ - l))vi{K,a,ii) = 0. (A.IO) 

We again expand the wavefunction over the Legendre polynomials 

vi{K,a,iJ,) = ^a£{K,a)Pi{iJ,). (A. 11) 

By iterating eq. (A. 4) twice, we obtain the following five-term recursion 

~ " [W^^mT^f'-' + (2^-1) (2^ + 3)"' + (2^ + 3)(2^ + 5)"^+V " 

(A.12) 

Eqs. (A. 6, 7, 9) still hold true. 

We finally consider the wave equation 

(A-aii)vi(a,ii,(p) = 0, (A.13) 

where A is the full Legendre operator, defined in eq. (4.9). Since the potential does not 
involve the azimuthal angle (p explicitly, we look for a solution vi proportional to e*"^*^, 
with m > being an integer. It is now natural to expand the function vi{a,iJ,,(fi) in the 
Legendre functions Pi^rn{l-i')j namely 

v^ia,fi,<f) = e^-^ (^£,MPi,M- (A.14) 

e>m 

These functions obey 

A(P,,^(//)e^-^) = -£{£+l)P,^^{i^)e^^^, (A.15) 
and the product //P^^^(/i) has the following expression 

(2£ + l)fiPe,M = (^ + 1 - m)Pe+i,M + {£ + m)Pe-i,M, (A.16) 
so that eq. (A.13) amounts to the three-term recursion relation 

(£ -\- fYl \ £ — Ul \ 

2£ + 3 Qw,m + ^a^-i.m j = 0. (A.17) 

The recursion equations (A. 5, 12, 17) are easily implemented numerically. 
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Appendix B. Semi-classical analysis 

The outcomes of the semi-classical analysis presented in this Appendix are used at 
various places in section 4. We first consider the Schrodinger equation (4.23). For large 
values of the complex parameter cr, we look for rapidly varying solutions of the form 

u ^ ^{xy^/"^ exp J $(y)dy, (B.l) 

with 

= aV{x). (B.2) 

This approach is analogous to the W.K.B. approximation in Quantum Mechanics, since 
the condition cr ^ 1 is equivalent to h being small. 

Because the potential V{x) is odd, we can restrict the analysis to the domain Re o" > 0. 
We introduce the notation 

K = y^. (B.3) 
Eq. (B.2) has real solutions for x > 0. We set 



p{x) = \/V{x) {x > 0), (B.4) 
so that ^{x) = Kp{x). Eq. (B.l) thus yields the basis of functions 

^±(^) ~ 7 ^TT72 ( ± KI{x)), (B.5) 

with 

l>oo 

I{x) = / p{y)dy {x > 0). (B.6) 

J X 

The above functions u± (x) are exponentially blowing up or decaying. The domains x > 
and cr > (and similarly a; < and cr < 0) are said to be classically forbidden. 

On the other hand, for x < 0, eq. (B.2) has imaginary solutions. We set 

q{x) = y/-V{x) {x < 0), (B.7) 
so that $(a;) = iKq{x). Eq. (B.l) thus yields the basis of functions 

u±(x) ^ — cxp ( ± iKI(x)), (B.8) 

{Kq{x)f^ V V 

with 

I{x)^ [ q{y)dy {x<0). (B.9) 
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The above functions u±{x) are oscillating and bounded, up to the prefactor in q(x)~^^'^. 
The domains a; < and u > (and similarly a; > and u < 0) are said to be classically 
allowed. 

The difficulty of the semi-classical analysis comes from the existence of three turning 
points, namely x — and x — > ±oo, where the momentum variable p{x) or q{x) vanishes. 
The estimates (B.5, 8) loose their meaning in the vicinity of the turning points, where a 
more careful analysis is required, to be presented now. 

We first investigate the basis of functions {wi,W2} for a; < 0. For x — oo and 
K — > oo, the Schrodinger equation (4.23) assumes the simpler form 



u" + {2Ke''Yu = 0. 



(B.IO) 



A basis of solutions to this equation is given by the Bessel functions Jq{z) and No{z), with 
z — 2Ke^ being a scaling variable. The boundary conditions (4.24) have to match the 
known small- 2; behavior of the Bessel functions, hence 



uiix) Ri Jo(2i^e^), 

U2ix) ^ (7r/2)Aro(2Ke^) -(\nK + 7^) Jo(2Ke^) (x 



-00), 



(B.ll) 



where 7^ denotes Euler's constant. The known large-^ behavior of the Bessel functions 
fixes the amplitudes of the integrals in eq. (B.8) for ui and U2, namely 



ui{x) 



nKq^x) 



1/2 



COS {KI{x) - 7r/4), 



/ 2 

U2{x) ^ [j^^j [(7r/2) sin {KI{x) - 7r/4) - (InK + 7^) cos {KI{x) - tt/A)] . 

(B.12) 

On the other hand, for a; — > and K — > 00, the Schrodinger equation (4.23) assumes the 
simpler form 

u" + K'^xu = {), (B.13) 

which is equivalent to Airy's equation. A basis of solutions is given by the Airy functions 
Ai(2;) and B\{z), with z = K'^^^x being again a scaling variable. The known behavior for 
2; — > — 00 of the Airy functions has to match eq. (B.12), hence 

ui(x) ^ 2^/2^-1/3 Jsin(K/) Ai{K^/^x) + cos(KI) Bi{K^/^x) 
U2{x) ^ 2^/2^-^/2 {(7^/2) - cos{KI) Ai{K^/^x) + sm{KI) Bi{K^/^x) (B.14) 



-(In K + 7s) \sm{KI) Ai{K'^/^x) + cos{KI) Bi{K'^/^ 



x 



for a; — > 0, with 
/ = /(O) = 



p(a;)da; 



AjJL 
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1/2 



¥ r(3/4) 

2 r(5/4) 



1.198140, 



(B.15) 



this definition being consistent with eqs. (B.6, 9). 

We now investigate in a similar way the functions {vi, ^2} for a; > 0. For a; — > 00 and 
— >^ 00, a basis of solutions is given by the modified Bessel functions Io{z) and Ko{z), with 
z = 2Ke~^. The boundary conditions (4.25) have to match the known small- 2; behavior 
of the Bessel functions, hence 

^ ^' B.16) 

V2Rii^o(2i^e-^) + (lnK + 7s)/o(2i^e-^) {x ^ 00). ^ ^ 

The known large-2 behavior of the Bessel functions only fixes the amplitude of the solution 
u+ of eq. (B.5) for vi and V2, namely 

V2{x) ^ {InK + ^e)vi{x) (a; > 0). 
Finally, the known behavior of the Airy functions as 2; ^ 00 yields 



vi{x)^2^/^K-y^Ai{K^/^x)e''', 
V2{x) P:: {InK + ^e)vi{x) (a; ^ 0). 



(B.18) 



The above expressions (B.14, 18) of both bases of solutions as a; — > and K ^ 00 
allow us to derive the following semi-classical estimates for the elements of the transfer 
matrix introduced in eq. (4.29) 

F{a) « [sin(i^/) - (2/7r)(lnK + 7^) cos(K/)]e-^^ 
Gia) ^ -(2/7r) cos(K7)e^^ 
H{a)^ {In K + ^^)F{a), 
F(-a) ^ (In K + ^^)G(a) (a = ^ 00). 

The estimate for G{a) directly yields the following expressions for its factors P(±(t), 
defined in eq. (4.41) 

KI 

P{a) ^ (3/7r)V2 

^ (B.20) 
Pi-a) ^ -2(3A)V2S^^ ia = K'^ 00). 

We also obtain from eq. (B.19) an estimate of the eigenvalues cr^ = K^, which are 
the zeroes of G{a), in the form 

{n + l/2)j (n>l). (B.21) 
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This semi-classical formula gives a very accurate description of the whole spectrum of the 
Schrodingcr equation (4.23). Indeed the relative error is maximal for the first non-zero 
eigenvalue, for which eq. (B.21) predicts Ki fti 3.933086, i.e., some 3.2% above the exact 
numerical value Ki = 3.811562. 

The semi- classical expression of the squared norms can be evaluated by inserting 
the estimates (B.19) into the identity (C.6). We thus obtain 

Nn ^ i-e''^-^ (n > 1). (B.22) 

In section 4 we also need the expression of the Mellin transform of the Airy function 
Ai(a;), namely 

POO 

m(s) = / x^Ai{x)dx (Res>-1), (B.23) 
Jo 

which we have not found in standard handbooks. The Airy equation implies the functional 
equation 

to(s + 3) = (s + l)(s + 2)m(s), (B.24) 
whose correctly normalised solution is 

We now consider the deformed Schrodingcr equation (4.63), where both a and k are 
non-zero. It turns out that only the spectrum of that wave equation will be needed. Hence 
we can content ourselves with the Sommerfeld quantisation formula. A similar treatment 
is used in section 4.5 for the full Legendre operator. 

The Sommerfeld formula reads 
j^^ ( - aV{x) - K^W{x)y^^dx = j^^ i^Y^ ~ ' ~ + V2) TT, (B.26) 

the integral being extended over the classically allowed domain, where the square-root is 
real. 

The implicit equation (B.26) for the semi-classical estimate of the eigenvalues ±a"„(/?) 
can be investigated in several limiting cases of interest. For small values of k, the integrand 
can be expanded in a straightforward way. We thus obtain 

TT / 2k^ \ 

(7n(«)^/' « Kn{n) = (n + 1/2) + ^^^^^^i^y + " " " j (n » 1, « « 1). (B.27) 

This expression confirms the general result (4.65). On the other hand, for cr = 0, it can be 
deduced from eq. (B.26) that the spheroidal equation (4.70) has imaginary eigenvalues of 
the form k = ±z^^, asymptotically given by the semi-classical estimate 

{n + l/2)n/2 (n > 1). (B.28) 
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Appendix C. Useful identities on the mixed Wronskians 

In this Appendix, we derive the identity (C.6) used in section 4, and more generally 
we give alternative expressions for the derivatives with respect to the spectral variable a 
of the mixed Wronskians -F((7), G((7), and H{a), which enter the transfer matrix (4.29). 

To do so, we start by considering the derivatives 

Ua{cr, X) = — {a = 1, 2), (C.l) 

which obey the inhomogeneous Schrodinger equation 

U^{a,x) - aV{x)Uc,{a,x) ^V{x)uc,{a,x) (a = 1,2). (C.2) 

These equations can be solved explicitly by varying the constants, along the lines of sections 
4.3 and 4.4. We thus obtain 

ui{a,y)u2{(T,y)V{y)dy + U2{(T,x) / ul{a,y)V{y)dy, 

-oo J — oo 



-oo 
rx 



/X px 
ul{a,y)V{y)dy + U2{(7,x) / ui{a,y)u2{(T,y)V{y)dy. 
-oo J — oo 



(C.3) 



By taking the a; — > oo limit of the above expressions, and using the asymptotic be- 
havior (4.30), we get the following expressions 



dF{a) 

da 
dG{a) 

da 
dH{a) 
da 
dF{-a) 
da 

with the definition 



G{a)N22{a) + F{a)N^2{<j), 
-F{a)Ni:,{a)-G{a)Ni2{a), 
-Fia)N-,^ia)-Gia)N-,2ia), 
-F{-a)Ni2{a) - H{a)Nii{a), 



(C.4) 



/oo 
Ua{a,x)ui3{a,x)V{x)dx (a,/? = 1,2). (C.5) 
-oo 

On the spectrum, i.e., for a = cr^, we have G{an) = 0, by definition. Furthermore eq. 
(4.34) implies Nii{an) = Nn/F'^{an), hence the identity 



- = -NnF[-an), (C.6) 



da J^^^^ F{a^) 

with Nn being the squared norm of the eigenfunction vi{anix), defined in eq. (4.35). 
The identity (C.6) is very general. It also holds for the K-dependent Schrodinger equation 
(4.63). 
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CAPTIONS FOR TABLES AND FIGURES 



Table 1: Conventions and notations for kinematic and other useful quantities. 

Table 2: Comparison of the numerical values of various quantities of interest from the 
exact solutions in the absence of internal reflections. First row: isotropic scattering, after 
ref. [15]; second row: very anisotropic scattering (this work). tqI* is the thickness of a 
skin layer; ti(1) and 7(1, 1) respectively yield the transmitted and reflected intensities in 
the normal direction; B{0) is the peak value of the enhancement factor at the top of the 
backscattering cone; t*AQ = ki£*AO is the dimensionless width of this backscattering 
cone. The third row gives the relative difference of the second case with respect to the 
flrst one. 

Figure 1: Laws of geometrical optics for the refraction of light by a large dielectric sphere. 

Figure 2: Plot of exact expressions for Ti{fi) in the absence of internal reflections. Dashed 
line: isotropic scattering, after ref. [15]. Full line: very anisotropic scattering (this work). 
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Table 1 





outside medium 


inside medium 


optical index 




77-0 = mni 


wavenumber 


ki — niLo/c — 2n/Xi 


ko — nooo/c — mki — 2n/Xo 


incidence angle 


e 


6' 


parallel wavevector 


p — ki cos 9 
= /co\/a*^ — 1 + 1/ 


P = ko cos e' 
= koi^ 


total reflection 
condition 


m < 1 and sin 6* > m 
(i.e. P imaginary) 


m > 1 and sin 6*' > 1/m 
(i.e. p imaginary) 



transverse 
wavevector 


|q| 


= q — ki sin 9 = ko sin 9' — ko^ 1 — jj? 


azimuthal angle 




reflection 
and 
transmission 
coefficients 


partial 
reflection 

total 
reflection 


^_ (P-pV _ ^///'-l + lA"^'^ 

^ 4Pp A|Ji^J -1 + 

\ R=l 
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Table 2 





To 




7(1,1) 


B{0) 


T*AQ 


isotropic 

(t* — ^) 


0.710446 


5.036475 


4.227681 


1.881732 


1/2 


very anisotropic 
(r* > 1) 


0.718211 


5.138580 


4.889703 


2 


0.555543 


A(%) 


1.1 


2.0 


15.7 


6.3 


11.1 
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